Thermodynamic Order Parameters and Statistical-Mechanical Measures for Characterization of the Burst and Spike Synchronizations of Bursting Neurons

Thermodynamic Order Parameters and Statistical-Mechanical Measures for Characterization of the Burst and Spike Synchronizations of Bursting Neurons

Abstract

We are interested in characterization of population synchronization of bursting neurons which exhibit both the slow bursting and the fast spiking timescales, in contrast to spiking neurons. Population synchronization may be well visualized in the raster plot of neural spikes which can be obtained in experiments. The instantaneous population firing rate (IPFR) , which may be directly obtained from the raster plot of spikes, is often used as a realistic collective quantity describing population behaviors in both the computational and the experimental neuroscience. For the case of spiking neurons, realistic thermodynamic order parameter and statistical-mechanical spiking measure, based on , were introduced in our recent work to make practical characterization of spike synchronization. Here, we separate the slow bursting and the fast spiking timescales via frequency filtering, and extend the thermodynamic order parameter and the statistical-mechanical measure to the case of bursting neurons. Consequently, it is shown in explicit examples that both the order parameters and the statistical-mechanical measures may be effectively used to characterize the burst and spike synchronizations of bursting neurons.

pacs:
87.19.lm, 87.19.lc

I Introduction

In recent years, brain rhythms which are observed in scalp electroencephalogram and local field potentials have attracted much attention (1). These brain rhythms emerge via synchronization between individual neuronal firings. Synchronization of firing activities may be used for efficient sensory and cognitive processing (e.g., feature integration, selective attention, and memory formation) (2); (3); (4). This kind of neural synchronization is also correlated with pathological rhythms associated with neural diseases such as epilepsy, Parkinson’s disease, and Alzheimer’s disease (5); (6); (7). Here, we are interested in characterization of these synchronous brain rhythms.

There are two basic types of neuronal firing activities, spiking and bursting (8). We are concerned about synchronization of bursting neurons. Bursting occurs when neuronal activity alternates, on a slow timescale, between a silent phase and an active (bursting) phase of fast repetitive spikings (9); (10); (11); (12); (13). In neural information transmission, burst input is more likely to have a stronger impact on the postsynaptic neuron than single spike input (14); (15); (16); (17). Intrinsically bursting neurons and chattering neurons in the cortex, thalamocortical relay neurons, thalamic reticular neurons, hippocampal pyramidal neurons, Purkinje cells in the cerebellum, pancreatic -cells, and respiratory neurons in pre-Botzinger complex are representative examples of bursting neurons (10); (13). These bursting neurons exhibit two different patterns of synchronization due to the slow and fast timescales of bursting activity. Burst synchronization (synchrony on the slow bursting timescale) refers to a temporal coherence between the active phase onset or offset times of bursting neurons, while spike synchronization (synchrony on the fast spike timescale) characterizes a temporal coherence between intraburst spikes fired by bursting neurons in their respective active phases (18); (19). Recently, many 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 burst and spike synchronizations) (20); (21); (22); (23); (24); (25); (26); (27); (28); (29); (30); (31); (32); (33); (34).

In this paper, we are concerned about practical characterization of the burst and spike synchronizations of bursting neurons. 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 collective quantity describing population behaviors in both the computational and the experimental neuroscience (2); (35); (36); (37); (38); (39); (40). We note that the experimentally-obtainable is in contrast to the ensemble-averaged potential which is often used as a population quantity in the computational and theoretical neuroscience, because to directly obtain in real experiments is very difficult. To overcome this difficulty, instead of , we employed as a population quantity, and developed realistic measures, based on , to make practical characterization of synchronization of spiking neurons in both the computational and the experimental neuroscience (41). The mean square deviation of plays the role of an order parameter used for characterizing synchronization transition of spiking neurons (42). The order parameter can be regarded as a “thermodynamic” measure because it concerns just the macroscopic quantity without considering any quantitative relation between and the microscopic individual spikes. Through calculation of , one can determine the threshold value for the spike synchronization. Moreover, to quantitatively measure the degree of spike synchronization, a “statistical-mechanical” spiking measure was introduced by taking into consideration both the occupation pattern and the pacing pattern of spikes in the raster plot. Particularly, the pacing degree between spikes was determined in a statistical-mechanical way by quantifying the average contribution of (microscopic) individual spikes to the (macroscopic) IPFR . Consequently, synchronization of spiking neurons may be well characterized in terms of these realistic thermodynamic order parameter and statistical-mechanical measure, and , based on .

The main purpose of our work is to characterize the burst and spike synchronizations of bursting neurons by extending the thermodynamic order parameter and the statistical-mechanical measure of spiking neurons (41) to the case of bursting neurons. For this aim, 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). Then, the mean square deviations of and play the role of realistic thermodynamic order parameters, and , used to determine the bursting and spiking thresholds for the burst and spike synchronization, respectively. We also consider another raster plot of bursting onset or offset times for more direct visualization of bursting behavior. From this type of raster plot, we can directly obtain the IPBR, or , without frequency filtering. Then, the time-averaged fluctuations of and also play the role of the order parameters, and , for the bursting transition. These bursting order parameters and are more direct ones than because they may be obtained directly without frequency filtering and they yield the same bursting threshold which is obtained through calculation of . As a next step, in the whole region of burst synchronization, the degree of burst synchronization seen in the raster plot of bursting onset or offset times may be well measured in terms of a statistical-mechanical bursting measure , introduced by considering both the occupation and the pacing patterns of bursting onset or offset times in the raster plot. In a similar way, we also develop a statistical-mechanical spiking measure , based on , to quantitatively measure the degree of the intraburst spike synchronization. Consequently, through separation of the slow bursting and the fast spiking timescales, burst synchronization may be well characterized in terms of both the bursting order parameters (, and ) and the statistical-mechanical bursting measure (), while characterization of intraburst spike synchronization can be made well by using the spiking order parameter () and the statistical-mechanical spiking measure (). To our knowledge, no measures characterizing intraburst spike synchronization of bursting neurons seem to be introduced previously. Hence, and are new realistic measures characterizing the intraburst spike synchronization. For the case of burst synchronization, our bursting order parameters (, and ) and the statistical-mechanical bursting measure () are also in contrast to the conventional measures such as the normalized order parameter (43); (44); (45); (46) and the burst phase order parameter (20); (23); (25); (47); (48). The normalized order parameter is given through dividing the order parameter (i.e., the time-averaged fluctuation of the ensemble-averaged potential ) by the average of time-averaged fluctuations of individual potentials. Our bursting order parameters, based on experimentally-obtainable IPBRs, are realistic ones when compared to , based on , because is very difficult to obtain in experiments. Furthermore, since shows both the bursting and spiking activities, plays the role of an order parameter for the “whole” synchronization (including both the burst and spike synchronizations) of bursting neurons, which is in contrast to our bursting order parameters characterizing just the burst synchronization. On the other hand, the burst phase order parameter is a “microscopic” measure quantifying the degree of coherence between (microscopic) individual burst phases without any explicit relation to the macroscopic occupation and pacing patterns of bursting onset or offset times visualized well in the raster plot, in contrast to our statistical-mechanical bursting measure . For our case of , the pacing degree (between the bursting onset or offset times) is determined in a statistical-mechanical way by taking into consideration the average of contributions of microscopic individual bursts to the macroscopic IPBR.

This paper is organized as follows. In Sec. II, as an example for characterization we describe an inhibitory population of bursting Hindmarsh-Rose (HR) neurons (49); (50); (51); (53). In Sec. III, through separation of the slow and fast timescales, we develop realistic order parameters and statistical-mechanical measures, based on IPBR and IPSR, which are applicable in both the computational and experimental neuroscience. Their usefulness for characterization of the burst and spike synchronizations is shown in explicit examples of bursting HR neurons. Finally, a summary is given in Section IV.

Ii Inhibitory Population of Bursting Hindmarsh-Rose Neurons

As an example for characterization, we consider an inhibitory population of globally-coupled bursting neurons. As an element in our coupled neural system, we choose the representative bursting HR neuron model which was originally introduced to describe the time evolution of the membrane potential for the pond snails (49); (50); (51); (53). The population dynamics in this neural network is governed by the following set of ordinary differential equations:

(1)
(2)
(3)
(4)

where

(5)
(6)

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

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 [Fig. 1(a)] to a bursting state [Fig. 1(b)]. As shown in Fig. 1(c), projection of the phase flow onto the plane seems to be a hedgehog-like attractor. Bursting activity [alternating between a silent phase and an active (bursting) phase of fast repetitive spikings] occurs on the hedgehog-like attractor [the body (spines) of the hedgehog-like attractor corresponds to the silent (active) phase]. Here, we consider the suprathreshold case of where each HR neuron exhibits spontaneous bursting activity without noise.

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) (45); (54). 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 (55). 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 ((56); (57)). Hence, rises fast and decays slowly.

Numerical integration of Eqs. (1)-(4) is done using the Heun method (58) (with the time step ms). For each realization of the stochastic process, we choose a random initial point for the th neuron with uniform probability in the range of , , , and .

Iii Characterization of The Burst and Spike Synchronizations in Terms of Thermodynamic Order Parameters and Statistical-Mechanical Measures

In this section, we extend the order parameter and the statistical-mechanical measure to the case of bursting neurons for characterization of population synchronization of bursting neurons. For this aim, we separate the slow and fast timescales of the bursting activity via 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 thermodynamic order parameters and statistical-mechanical measures, based on and , and show their usefulness for characterization of the burst and spike synchronizations in explicit examples of HR bursting neurons.

As an example for characterization, we consider an inhibitory population of globally-coupled bursting HR neurons for and set the coupling strength as . By varying the noise intensity , we characterize the population synchronization of bursting HR neurons. In computational neuroscience, an ensemble-averaged global potential,

(7)

is often used for describing emergence of population synchronization. Throughout this study, we consider the population behaviors after the transient time of ms. Figure 2(a) shows an oscillating global potential for a synchronous case of . For comparison, an individual potential of the 1st HR neuron is also shown in Fig. 2(a). In contrast to , each HR neuron fires sparse burstings about once per three global cycles of . An active phase of the bursting activity begins at a bursting onset time and ends at a bursting offset time. At the bursting onset (offset) time, each individual potential of the th bursting neuron passes the threshold of from below (above); the bursting onset (offset) times of the 1st neuron (after the transient time) are denoted by the solid (open) circles in Fig. 2(a). Although the global potential is an important population-averaged quantity to describe synchronization in computational neuroscience, it is practically difficult to directly obtain 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); (35); (36); (37); (38); (39); (40). The IPFR is obtained from the raster plot of neural spikes which is a collection of spike trains of individual neurons. Such raster plots of spikes, where population spike synchronization may be well visualized, are fundamental data in experimental neuroscience. As examples of population bursting states, Figs. 2(b1)-2(b5) show the raster plots of neural spikes for various values of noise intensity : synchronized bursting states for 0.01, 0.04, and 0.06, and unsynchronized bursting state for . To obtain a smooth IPFR from the raster plot of spikes, we employ the kernel density estimation (kernel smoother) (59). Each spike in the raster plot is convoluted (or blurred) with a kernel function to obtain 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(c1)-2(c5) show smooth IPFR kernel estimates of band width ms. For , clear “bursting bands,” each of which is composed of “stripes” of spikes, appear successively at nearly regular time intervals [see Fig. 2(b1)]; a magnification of the 1st bursting band is given in Fig. 6(a1). For this case of , in addition to burst synchronization [synchrony on the slow bursting timescale ( ms)], spike synchronization [synchrony on the fast spike timescale ms)] occurs in each bursting band. As a result of this complete synchronization, the IPFR kernel estimate exhibits a bursting activity [i.e., fast spikes appear on a slow wave in ], as shown in Fig. 2(c1). However, as is increased, loss of spike synchronization begins to occur in each bursting band due to a smearing of spiking stripes. As an example, see the case of where the raster plot of spikes and the IPFR kernel estimate are given in Figs. 2(b2) and 2(c2), respectively. Smearing of the spiking stripes is well seen in the magnified 1st bursting band of Fig. 6(a3). Hence, the amplitude of decreases, as shown in Fig. 2(c2). As is further increased and passes a spiking noise threshold , complete loss of spike synchronization occurs in each bursting band (i.e., spikes become incoherent within each bursting band). Consequently, only the burst synchronization (without spike synchronization) occurs [e.g., see the case of in Figs. 2(b3) and 2(c3)]. For this case, shows a slow-wave oscillation without spikes. With increasing , such “incoherent” bursting bands become more and more smeared, and thus the degree of burst synchronization decreases [e.g., see Figs. 2(b4) and 2(c4) for ]. Consequently, the amplitude of is further decreased. With further increase in , incoherent bursting bands begin to overlap, which eventually results in the complete loss of burst synchronization when passing a larger bursting noise threshold . In this way, completely unsynchronized states with nearly stationary appear, as shown in Figs. 2(b5) and 2(c5) for .

We note that the IPFR kernel estimate is a population quantity describing the “whole” combined collective behaviors (including both the burst and spike synchronizations) of bursting neurons. For more clear investigation of population synchronization, we separate the slow bursting timescale and the fast spiking timescale via frequency filtering, and decompose the IPFR kernel estimate into the IPBR and the IPSR . Through low-pass filtering of with cut-off frequency of 10 Hz, we obtain the regularly-oscillating IPBR (containing only the bursting behavior without spiking) in Figs. 2(d1)-2(d5). As is increased, the amplitude of decreases gradually, and eventually becomes nearly stationary when passing the bursting noise threshold . For more direct visualization of bursting behavior, we consider another raster plot of bursting onset or offset times, from which we can directly obtain the IPBR kernel estimate of band width ms, or , without frequency filtering. Figures 2(e1)-2(e5) show the raster plots of the bursting onset times, while the raster plots of the bursting offset times are shown in Figs. 2(f1)-2(f5). From these raster plots of the bursting onset (offset) times, we obtain smooth IPBR kernel estimates, [] in Figs. 2(g1)[(h1)]-2(g5)[(h5)]. For , clear “bursting stripes” [composed of bursting onset (offset) times and indicating burst synchronization] are formed in these raster plots; the bursting onset and offset stripes are time-shifted [see Figs. 2(e1) and 2(f1)]. The corresponding IPBR kernel estimates, and , for show regular oscillations with the same population bursting frequency Hz), as shown in Figs. 2(g1) and 2(h1), although they are phase-shifted. As is increased, the bursting stripes in the raster plots become smeared and begin to overlap, and thus the degree of the burst synchronization decreases. Consequently, the amplitudes of both and decrease gradually (e.g., see the cases of , 0.04, and 0.06). Eventually, when passing 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. 2(g5) and 2(h5) for . In this way, and are more direct ones for describing the bursting behaviors than .

As is well known, a conventional order parameter, based on the ensemble-averaged global potential , is often used for describing transition from asynchrony to synchrony in computational neuroscience (43); (44); (45); (46). Here, instead of which is very difficult to obtain in experiments, we use an experimentally-obtainable IPBR (which is obtained from the IPFR kernel estimate via low-pass filtering), and develop a realistic bursting order parameter for the bursting transition, which may be applicable in both the computational and the experimental neuroscience. The mean square deviation of ,

(10)

plays the role of a bursting order parameter , where the overbar represents the time average. The order parameter may be regarded as a thermodynamic measure because it concerns just the macroscopic IPBR without any consideration between and microscopic individual burstings. Here, we discard the first time steps of a trajectory as transients for ms, and then we compute by following the trajectory for ms. As is increased, exhibits more regular oscillations for the case of burst synchronization, while becomes more stationary for the case of burst unsynchronization. Hence, the bursting order parameter , representing time-averaged fluctuations of , approaches a non-zero (zero) limit value for the synchronized (unsynchronized) bursting state in the thermodynamic limit of . Figure 3(a) shows plots of the order parameter versus . For ), synchronized bursting states exist because the values of become saturated to non-zero limit values. As passes the bursting noise threshold , the bursting order parameter tends to zero as , and hence a transition to unsynchronized bursting states occurs because the noise spoils the burst synchronization. In addition to , we also employ another IPBR kernel estimates, and , (which are directly obtained from the raster plots of bursting onset and offset times), and introduce realistic bursting order parameters, and :

(11)

Figures 3(b) and 3(c) show plots of the bursting order parameters and versus , respectively. Like the case of , when passing the same bursting noise threshold , the bursting order parameters and also go to zero as , and hence a transition to burst unsynchronization occurs. In this way, the noise threshold for the bursting transition may be well determined through calculation of the three realistic bursting order parameters, , and . Particularly, and are more direct ones than because they are based on the IPBRs and which are directly obtained from the raster plots of the bursting onset and offset times without frequency filtering, respectively.

As a next step, within the burst-synchronized region (), we measure the degree of burst synchronization in terms of a realistic statistical-mechanical bursting measure , based on the IPBR kernel estimates and . Previously, a statistical-mechanical spiking measure, based on the ensemble-averaged global potential , was developed for characterization of spike synchronization of spiking neurons (60). However, the spiking measure, based on , is practically inapplicable to the case of experimental neuroscience because to obtain in experiments is difficult. Hence, instead of , we used the experimentally-obtainable IPSR kernel estimate, and developed a refined version of statistical-mechanical spiking measure, based on the IPSR, to characterize spike synchronization of spiking neurons in both the experimental and the computational neuroscience (41). Here, we extend the realistic statistical-mechanical measure of spiking neurons (based on the IPSR) to the case of bursting neurons for measurement of the degree of the burst synchronization. As shown in Figs. 2(e1)[(f1)]-2(e5)[(f5)], burst synchronization may be well visualized in the raster plots of bursting onset (offset) times. For the synchronous bursting case, bursting stripes (composed of bursting onset (offset) times and indicating population burst synchronization) appear in the raster plots. As an example, we consider a synchronous bursting case of . The raster plot in Fig. 4(a1) is composed of partially-occupied and smeared stripes of bursting onset times. To measure the degree of burst synchronization seen in the raster plot, we develop a statistical-mechanical bursting onset measure , based on , by considering both the occupation pattern and the pacing pattern of the bursting onset times in the bursting onset stripes. The bursting onset measure of the th bursting onset stripe is defined by the product of the occupation degree of bursting onset times (representing the density of the th bursting onset stripe) and the pacing degree of bursting onset times (denoting the smearing of the th bursting onset stripe):

(12)

The occupation degree of bursting onset times in the th bursting onset stripe is given by the fraction of HR neurons which fire burstings:

(13)

where is the number of HR neurons which fire burstings in the th bursting onset stripe. For the full occupation , while for the partial occupation . The pacing degree of bursting onset times in the th bursting onset stripe can be determined in a statistical-mechanical way by taking into account their contributions to the macroscopic IPBR kernel estimate . Figure 4(a2) shows a time series of the IPBR kernel estimate for ; local maxima and minima are denoted by solid and open circles, respectively. Obviously, central maxima of between neighboring left and right minima of coincide with centers of bursting onset stripes in the raster plot. The global bursting onset cycle starting from the left minimum of which appears first after the transient time ms) is regarded as the 1st one, which is denoted by . The 2nd global bursting onset cycle begins from the next following right minimum of , and so on. (The 1st global bursting onset cycle corresponds to the 2nd bursting onset stripe in Fig. 2(e1) because the minimum of the global bursting onset cycle, corresponding to the 1st bursting onset stripe in Fig. 2(e1), lies for ms.) Then, we introduce an instantaneous global bursting onset phase of via linear interpolation in the two successive subregions forming a global bursting onset cycle (41); (60); (61), as shown in Fig. 4(a3). The global bursting onset phase between the left minimum (corresponding to the beginning point of the th global bursting onset cycle) and the central maximum is given by:

(14)

and between the central maximum and the right minimum (corresponding to the beginning point of the th global bursting onset cycle) is given by

(15)

where is the beginning time of the th global bursting onset cycle (i.e., the time at which the left minimum of appears in the th global bursting onset cycle) and is the time at which the maximum of appears in the th global bursting onset cycle. Then, the contribution of the th microscopic bursting onset time in the th bursting onset stripe occurring at the time to is given by , where is the global bursting onset phase at the th bursting onset time [i.e., ]. A microscopic bursting onset time makes the most constructive (in-phase) contribution to when the corresponding global onset phase is (), while it makes the most destructive (anti-phase) contribution to when is . By averaging the contributions of all microscopic bursting onset times in the th bursting onset stripe to , we obtain the pacing degree of bursting onset times in the th bursting onset stripe:

(16)

where is the total number of microscopic bursting onset times in the th bursting onset stripe. By averaging of Eq. (12) over a sufficiently large number of bursting onset stripes, we obtain the realistic statistical-mechanical bursting onset measure , based on the IPSR kernel estimate :

(17)

For we follow bursting onset stripes and get , , and in each th bursting onset stripe, which are shown in Figs. 4(c1), 4(d1), and 4(e1), respectively. Due to sparse burstings of individual HR neurons, the average occupation degree (=, where denotes the average over bursting onset stripes, is small. Hence, only a fraction (about 1/3) of the total HR neurons fire burstings in each bursting onset stripe. On the other hand, the average pacing degree (= is large in contrast to . Consequently, the realistic statistical-mechanical bursting onset measure (=), representing the degree of burst synchronization seen in the raster plot of bursting onset times, is about 0.31. The main reason for the low degree of burst synchronization is mainly due to partial occupation. In this way, the realistic statistical-mechanical bursting onset measure can be used effectively to measure the degree of burst synchronization because concerns not only the pacing degree, but also the occupation degree of bursting onset times in the bursting onset stripes of the raster plot.

In addition to the above case of bursting onset times, we also consider population burst synchronization between the bursting offset times. Figures 4(b1) and 4(b2) show the raster plot composed of two stripes of bursting offset times and the corresponding IPBR for , respectively; the 1st and 2nd global bursting offset cycles, and , are shown. [Since the 1st global cycle of offset times, , follows the 1st global cycle of onset times, , the 1st bursting offset stripe in Fig. 4(b1) corresponds to the 2nd bursting offset stripe in Fig. 2(f1).] Then, as in the case of , one can introduce an instantaneous global bursting offset phase of via linear interpolation in the two successive subregions forming a global bursting offset cycle [see Fig. 4(b3)]. Similar to the case of bursting onset times, we measure the degree of burst synchronization seen in the raster plot of bursting offset times in terms of a statistical-mechanical bursting offset measure , based on , by considering both the occupation pattern and the pacing pattern of the bursting offset times in the bursting offset stripes. The bursting offset measure in the th bursting offset stripe also is defined by the product of the occupation degree of bursting offset times and the pacing degree of bursting offset times in the th bursting offset stripe. We also follow bursting offset stripes and get , , and in each th bursting offset stripe for , which are shown in Figs. 4(c2), 4(d2), and 4(e2), respectively. For this case of bursting offset times, (=, (=, and (=) . The pacing degree of offset times is a little smaller than the pacing degree of the onset times (), although the occupation degrees of the onset and offset times are the same. We take into consideration both cases of the onset and offset times equally and define the average occupation degree , the average pacing degree , and the statistical-mechanical bursting measure as follows:

(18)

By increasing the noise intensity , we follow 500 bursting onset and 500 bursting offset stripes and characterize burst synchronization in terms of (average occupation degree), (average pacing degree), and (statistical-mechanical bursting measure) for 15 values of in the whole region of burst synchronization ], and the results are shown in Figs. 5(a)-5(c). As is increased, the average occupation degree (denoting the average density of bursting stripes in the raster plot) decreases very slowly around , because a little tendency for noise-induced skipping of burstings in individual HR neurons occurs (45). On the other hand, with increasing , the average pacing degree (representing the average smearing of the bursting stripes in the raster plot) decreases rapidly due to destructive role of noise spoiling burst synchronization. The statistical-mechanical bursting measure also makes a rapid decrease because of a rapid drop in . Both and show quadratic decreases because they are well fitted with quadratic functions: and . In this way, we measure the degree of burst synchronization in terms of the realistic statistical-mechanical bursting measure in the whole synchronized region, and find that reflects the degree of burst synchronization seen in the raster plot of onset and offset times very well.

Unlike the case of spiking neurons (showing only the spike synchronization), bursting neurons may exhibit both the burst and the spike synchronizations. From now on, we investigate the intraburst spike synchronization of bursting HR neurons by varying the noise intensity . Figures 6(a1)-6(a6) and Figures 6(b1)-6(b6) show the raster plots of intraburst spikes and the corresponding IPFR kernel estimates during the 1st global bursting cycle of the low-pass filtered IPBR , respectively for various values of : synchronized spiking states for , 0.005, 0.01, and 0.02, and unsynchronized spiking states for and 0.08. As mentioned above, exhibits the whole combined population behaviors including the burst and spike synchronizations with both the slow bursting and the fast spiking timescales. Hence, through band-pass filtering of [with the lower and the higher cut-off frequencies of 30 Hz (high-pass filter) and 90 Hz (low-pass filer)], we obtain the IPSRs , which are shown in Figs. 6(c1)-6(c6). Then, the intraburst spike synchronization may be well described in terms of the IPSR . For , clear 8 “spiking stripes” (composed of spikes and indicating population spike synchronization) appear in the bursting band of the 1st global bursting cycle of the IPBR [see Fig. 6(a1)], and the IPFR kernel estimate exhibits a bursting activity [i.e., fast spikes appear on a slow wave in ] due to the complete synchronization (including both the burst and spike synchronizations), as shown in Fig. 6(b1). However, the band-pass filtered IPSR shows only the fast spiking oscillations (without a slow wave) with the population spiking frequency Hz) [see Fig. 6(c1)]. As is increased, spiking stripes in the bursting band become more and more smeared (e.g., see the cases of , 0.01, and 0.02). As a result, the amplitude of the IPSR decreases due to loss of spike synchronization. Eventually, when passing the spiking noise threshold , spikes become completely scattered within the bursting band (i.e., intraburst spikes become completely incoherent), and hence complete loss of spike synchronization occurs in the bursting band. As an example, see the case of . For this case, the IPSR becomes nearly stationary, while the IPFR kernel estimate shows a slow-wave oscillation (without spikes) due to the burst synchronization. Thus, for only the burst synchronization may occur. With further increase in , the incoherent bursting band expands, fills the whole region of the global bursting cycle, and overlaps with nearest bursting bands. Consequently, complete loss of burst synchronization also occurs when passing the larger bursting noise threshold . Thus, for completely unsynchronized states with nearly stationary appear (e.g., see the case of ).

For characterization of the intraburst spiking transition, we employ the experimentally-obtainable IPSR (which is obtained from the IPFR kernel estimate via band-pass filtering), and develop a realistic spiking order parameter , which may be applicable in both the computational and the experimental neuroscience. The mean square deviation of in the th global bursting cycle,

(19)

plays the role of a spiking order parameter in the th global bursting cycle of the low-pass filtered IPBR . By averaging over a sufficiently large number of global bursting cycles, we obtain the realistic spiking order parameter:

(20)

By following bursting cycles, we obtain the spiking order parameter . Figure 6(d) shows plots of versus . For ), synchronized spiking states exist because the values of become saturated to non-zero limit values in the thermodynamic limit of . However, when passing the spiking noise threshold , the spiking order parameter tends to zero as , and hence a transition to unsynchronized spiking states occurs because the noise spoils the spike synchronization. In this way, the spiking noise threshold may be well determined through calculation of the realistic spiking order parameter .

Finally, within the whole region of the intraburst spike synchronization (), we measure the degree of intraburst spike synchronization in terms of a realistic statistical-mechanical spiking measure , based on the IPSR . As shown in Figs. 6(a1)-6(a6), spike synchronization may be well visualized in the raster plot of spikes. For the synchronous spiking case, spiking stripes (composed of spikes and indicating population spike synchronization) appear within the bursting bands of the raster plot. As an example, we consider a synchronous spiking case of . Figures 7(a1) and 7(a2) show a magnified raster plot of neural spikes and the IPSR , corresponding to the 1st global bursting cycle of the low-pass filtered IPBR [bounded by a vertical dash-dotted lines: ]. Within the 1st global cycle, the bursting band [bounded by the vertical dotted lines: ], corresponding to the 1st global active phase, is composed of 8 stripes of spikes, as shown in Fig. 7(a1); (maximum of within the 1st global bursting cycle) is the global active phase onset time, and (maximum of within the 1st global bursting cycle) is the global active phase offset time. In the bursting band, the maxima (minima) of are denoted by solid (open) circles, and 8 global spiking cycles exist in the 1st global bursting cycle of , as shown in Fig. 7(a2). For , each th global spiking cycle , containing the th maximum of , begins at the left nearest-neighboring minimum of and ends at the right nearest-neighboring minimum of , while for both extreme cases of and 8, begins at [the beginning time of the 1st global bursting cycle of ] and ends at [the ending time of the 1st global bursting cycle of ]. Then, as in the case of the global bursting phase, one can introduce an instantaneous global spiking phase of via linear interpolation in the two successive subregions (the left subregion joining the left beginning point and the central maximum and the right subregion joining the central maximum and the right ending point) forming a global spiking cycle [see Fig. 7(a3)]. Similar to the case of burst synchronization, we measure the degree of the intraburst spike synchronization seen in the raster plot in terms of a statistical-mechanical spiking measure, based on , by considering both the occupation pattern and the pacing pattern of spikes in the global spiking cycles. The spiking measure of the th global spiking cycle in the 1st global bursting cycle is defined by the product of the occupation degree of spikes (representing the density of spikes in the th global spiking cycle) and the pacing degree of spikes (denoting the smearing of spikes in the th global spiking cycle). Figures 7(b1)-7(b3) show the plots of , , and , respectively. For the 1st global bursting cycle, the spiking-averaged occupation degree (=) , the spiking-averaged pacing degree (=) , and the spiking-averaged statistical-mechanical spiking measure (=) , where represents the average over the spiking cycles. We also follow bursting cycles and get , , and in each th global bursting cycle for , which are shown in Figs. 7(c1), 7(c2), and 7(c3), respectively. Then, through average over all bursting cycles, we obtain the bursting-averaged occupation degree (=, the bursting-averaged pacing degree (=, and the bursting-averaged statistical-mechanical spiking measure (= for . We note that , , and are obtained through double-averaging over the spiking and bursting cycles. When compared with the bursting case of and for , a fraction (about 3/4) of the HR neurons exhibiting the bursting active phases fire spikings in the spiking cycles, and the pacing degree of spikes is about 60 percentage of the pacing degree of burstings . Consequently, the statistical-mechanical spiking measure becomes about 45 percentage of the statistical-mechanical bursting measure for . We increase the noise intensity and obtain , , and . However, as will be seen below, with increasing , decreases very rapidly in an exponential way, in contrast to the bursting case. Hence, for more accurate results, we repeat the process to get , , and for multiple realizations. Thus, we obtain (average occupation degree of spikes in the global spiking cycles), (average pacing degree of spikes in the global spiking cycles), and (average statistical-mechanical spiking measure in the global spiking cycles) through average over all realizations. For each realization, we follow 100 bursting cycles, and obtain , , and via average over 20 realizations. Through these multiple-realization simulations, we characterize intraburst spike synchronization in terms of , , and for 8 values of in the whole region of spike synchronization [], which are shown in Figs. 7(d1)-7(d3), respectively. As is increased, the average occupation degree decreases very slowly around due to a little tendency for noise-induced subtracting of spikes in individual HR neurons, while the average pacing degree decreases very rapidly due to destructive role of noise spoiling spike synchronization. The average statistical-mechanical spiking measure also makes a rapid decrease because of a rapid drop in . Both and exhibit exponential decreases because they are well fitted with exponential functions: and . In this way, we measure the degree of intraburst spike synchronization in terms of the realistic statistical-mechanical spiking measure in the whole synchronized region, and find that reflects the degree of intraburst spike synchronization seen in the raster plot very well. Finally, we note that the exponential loss in the degree of spike synchronization is much faster than the quadratic loss in the degree of the burst synchronization. As a result, the break-up of the spike synchronization occurs first at the smaller spiking noise threshold , and then the burst synchronization disappears at the larger bursting noise threshold .

Iv Summary

We have extended the order parameter and the statistical-mechanical measure to the case of bursting neurons. Their usefulness for characterization of the burst and spike synchronizations has been shown in explicit examples of bursting HR neurons by varying the noise intensity . We note that population synchronization may be well visualized in the raster plot of neural spikes which may be obtained in experiments. Unlike the case of spiking neurons, bursting neurons show firing patterns with two timescales: a fast spiking timescale and a slow bursting timescale that modulates the spiking activity. Hence, the IPFR kernel estimate , which is obtained from the raster plot of spikes, shows collective behaviors with both the slow bursting and the fast spiking timescales. For our purpose, we separate the slow bursting and the fast spiking timescales via frequency filtering, and decompose the IPFR kernel estimate into the IPBR and the IPSR . Based on and , we have developed the bursting and spiking order parameters and which may be used to determine the bursting and spiking noise thresholds, and , for the burst and spike synchronizations. When passing and , loss of the burst and spike synchronizations occurs due to destructive role of noise spoiling the burst and spike synchronizations, respectively. As a next step, the degree of burst synchronization seen in the raster plot of bursting onset or offset times has been well measured in the whole region of burst synchronization in terms of a statistical-mechanical bursting measure , introduced by considering both the occupation and the pacing patterns of bursting onset or offset times in the raster plot. Similarly, we have also developed a statistical-mechanical spiking measure , based on , and measured the degree of the intraburst spike synchronization well. Thus, the statistical-mechanical bursting and spiking measures have been found to reflect both the occupation and the pacing degrees of bursting onset or offset times and spikes seen in the raster plot very well. Furthermore, it has also been found that the exponential loss in the degree of spike synchronization is much faster than the quadratic loss in the degree of the burst synchronization. Hence, the intraburst spike synchronization breaks up first at the smaller spiking noise threshold , and then the burst synchronization disappears at the larger bursting noise threshold . Consequently, we have shown in explicit examples that the order parameters and the statistical-mechanical measures may be effectively used to determine the bursting and spiking thresholds for the burst and the spike synchronizations and also to quantitatively measure the degree of the burst and the spike synchronizations, respectively.

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).
Figure 1: Single bursting HR neuron. Time series of the fast membrane potential for (a) and (b) . (c) Projection of the phase flow onto the plane for .
Figure 2: Population bursting states for various values of in an inhibitory ensemble of globally-coupled bursting HR neurons for and : synchronized bursting states for 0.01, 0.04, and 0.06, and unsynchronized bursting state for . (a) Time series of the ensemble-averaged global potential and time series of the individual potential of the 1st neuron for . The dotted horizontal line () represents the bursting threshold (the solid and open circles denote the bursting onset and offset times, respectively), while the dashed horizontal line () represents the spiking threshold within the active phase. Raster plots of neural spikes for (b1)-(b5), time series of IPFR kernel estimate for (c1)-(c5), time series of low-pass filtered IPBR (cut-off frequency=10 Hz) for (d1)-(d5), raster plots of active phase (bursting) onset times for (e1)-(e5), raster plots of active phase (bursting) offset times for (f1)-(f5), time series of IPBR kernel estimate for (g1)-(g5), and time series of IPBR kernel estimate for (h1)-(h5). The band width of the Gaussian kernel function is 1 ms for the IPFR kernel estimate and 50 ms for the IPBR kernel estimates and .
Figure 3: Determination of the bursting noise threshold for the bursting transition in terms of realistic thermodynamic bursting order parameters in an inhibitory ensemble of globally-coupled bursting HR neurons for and . Plots of bursting order parameters (a) [based on ], (b) [based on ], and (c) [based on ] versus . For each , we compute the bursting order parameters, , , and by following a trajectory for ms after discard the transients for ms.
Figure 4: Realistic statistical-mechanical bursting measures for measurement of the degree of burst synchronization, based on the IPBR kernel estimates, and , in an inhibitory population of ) globally-coupled bursting HR neurons for and in the case of . (a1) Raster plot of bursting onset times, (a2) time series of the IPBR kernel estimate , and (a3) the global bursting onset phase . (b1) Raster plot of bursting offset times, (b2) time series of the IPBR kernel estimate , and (b3) the global bursting offset phase . In (a2)-(a3) and (b2)-(b3), vertical dashed and dotted lines represent the times at which local minima and maxima (denoted by open and solid circles) of and