Accepted Manuscript Estimating the correlation between bursty spike trains and local field potentials Zhaohui Li, Gaoxiang Ouyang, Li Yao, Xiaoli Li PII: DOI: Reference:

S0893-6080(14)00110-5 http://dx.doi.org/10.1016/j.neunet.2014.05.011 NN 3334

To appear in:

Neural Networks

Received date: 30 May 2013 Revised date: 17 May 2014 Accepted date: 23 May 2014 Please cite this article as: Li, Z., Ouyang, G., Yao, L., & Li, X. Estimating the correlation between bursty spike trains and local field potentials. Neural Networks (2014), http://dx.doi.org/10.1016/j.neunet.2014.05.011 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

*Title Page (With all author details listed)

Estimating the correlation between bursty spike trains and local field potentials Zhaohui Lia, Gaoxiang Ouyangb,c, Li Yaob,c, Xiaoli Lib,c* aSchool

of Information Science and Engineering, Yanshan University, Qinhuangdao 066004, People’s

Republic of China bState

Key Laboratory of Cognitive Neuroscience and Learning and IDG/McGovern Institute for Brain

Research, Beijing Normal University, Beijing, China cCenter

for Collaboration and Innovation in Brain and Learning Sciences, Beijing Normal University,

Beijing, China

__________ * Corresponding author at: State Key Laboratory of Cognitive Neuroscience and Learning, Beijing Normal University, No.19, XinJieKouWai St., HaiDian District, Beijing, 100875, P. R. China E-mail address: [email protected] (Xiaoli Li). Tel. 86-010-5880 2032

*Manuscript Click here to view linked References

ABSTRACT

To further understand rhythmic neuronal synchronization, an increasingly useful method is to determine the relationship between the spiking activity of individual neurons and the local field potentials (LFPs) of neural ensembles. Spike field coherence (SFC) is a widely used method for measuring the synchronization between spike trains and LFPs. However, due to the strong dependency of SFC on the burst index, it is not suitable for analyzing the relationship between bursty spike trains and LFPs, particularly in high frequency bands. To address this issue, we developed a method called weighted spike field correlation (WSFC), which uses the first spike in each burst multiple times to estimate the relationship. In the calculation, the number of times that the first spike is used is equal to the spike count per burst. The performance of this method was demonstrated using simulated bursty spike trains and LFPs, which comprised sinusoids with different frequencies, amplitudes, and phases. This method was also used to estimate the correlation between pyramidal cells in the hippocampus and gamma oscillations in rats performing behaviors. Analyses using simulated and real data demonstrated that the WSFC method is a promising measure for estimating the correlation between bursty spike trains and high frequency LFPs.

Keywords: burst; correlation; local field potential; phase locking; spike train

1. Introduction

The advent of multi-electrode arrays has facilitated the simultaneous recording of the spiking activity of multiple neurons and neural ensembles, which provides an important method for investigating fundamental issues related to neural coding (Claverol-Tinture, Cabestany, & Rosell, 2007; Galashan, et al., 2011; Schwartz, 2004; Stafford, Sher, Litke, & Feldheim, 2009). The voltage signals obtained are generally separated into two types: the spikes or action potentials, which are fired by neurons and identified by high-pass filtering, detection, and sorting; and the local field potentials (LFPs), which are the total synaptic currents in the neuronal circuit and are obtained by low-pass filtering the original wideband signal (Mizuseki, Sirota, Pastalkova, & Buzsaki, 2009; Perelman & Ginosar, 2007). The interactions between the spikes of single neurons, i.e., spike trains, and the ongoing LFP oscillations are becoming hot topics in neuroscience because they allows us to study how the activities of individual neurons are related to those of the larger-scale networks in which they are embedded. Their significance has been shown to be associated with high-level brain functions, such as attention (Chalk, et al., 2010; Fries, Reynolds, Rorie, & Desimone, 2001), memory (Harris, et al., 2002; Le Van Quyen, et al., 2008; Lee, Simpson, Logothetis, & Rainer, 2005), motor tasks (Courtemanche, Pellerin, & Lamarre, 2002; Hagan, Dean, & Pesaran, 2012; van Wingerden, Vinck, Lankelma, & Pennartz, 2010), and sensory processing (Eggermont & Smith, 1995; Fries, Roelfsema, Engel, Konig, & Singer, 1997; Pienkowski & Eggermont, 2011; Xu, Jiang, Poo, & Dan, 2012).

To quantify the correlation between spike trains and LFPs, a wide variety of spike-LFP measures have been introduced in the past few years, e.g., the phase histogram, which is calculated by summing spikes that occur at different LFP phases (Csicsvari, Jamieson, Wise, & Buzsaki, 2003); the pairwise phase consistency, which is a bias-free measure of rhythmic neuronal synchronization (Vinck, Battaglia, Womelsdorf, & Pennartz, 2012; Vinck, van Wingerden, Womelsdorf, Fries, & Pennartz, 2010); phase locking, which is evaluated by applying the Rayleigh test for circular uniformity to the spike phase distribution (Colgin, et al., 2009; Siapas, Lubenov, & Wilson, 2005; Sirota, et al., 2008); and coherency, which is obtained by normalizing the cross-spectrum of each process (spike train and LFP) with the

spectrum of each process (Chalk, et al., 2010; Gregoriou, Gotts, Zhou, & Desimone, 2009; Jarvis & Mitra, 2001; Pesaran, Pezaris, Sahani, Mitra, & Andersen, 2002). In addition, another commonly used method for studying spike-LFP interactions is the spike field coherence (SFC), which measures the synchronization between spike trains and LFPs as a function of the frequency, where it takes values between 0% (complete lack of synchronization) and 100% (complete synchronization) (Fries, et al., 2001; Fries, et al., 1997). The SFC can be used to describe the strength of synchronization between spike times and a particular phase of the LFP oscillation at a certain frequency. The SFC has been employed to investigate memory formation in humans (Rutishauser, Ross, Mamelak, & Schuman, 2010), the neural mechanism of visual attention in macaque monkeys (Chalk, et al., 2010; Fries, et al., 2001), stimulus-specific synchronization in the primary visual cortex of awake cats exhibiting behaviors (Siegel & Konig, 2003), and other brain functions (Fries, Schroder, Roelfsema, Singer, & Engel, 2002; Issa & Wang, 2011; Lewandowski & Schmidt, 2011; Tiesinga, Fellous, Salinas, Jose, & Sejnowski, 2004; Wang, Iliescu, Ma, Josic, & Dragoi, 2011). However, this method also has the drawback that the estimated coherence is biased when only a small number of spikes are available (Grasse & Moxon, 2010; Vinck, Lima, et al., 2010; Vinck, van Wingerden, et al., 2010). Recently, a new pairwise measure that is not biased by the number of spikes was developed to address this problem (Vinck, et al., 2012). In addition to the spike number, the temporal structure of spike trains has an effect on the SFC calculation. Bursts, i.e., short episodes of high-frequency spike firing (Mammone & Morabito, 2008), are commonly observed structures in spike trains. In the present study, we modified the SFC method for bursty spike trains and LFPs to estimate correlations. Given that the coherence is generally used to measure the linear association between two signals in the frequency domain, we use the correlation to denote the relationships between spike trains and LFPs in the present study, which emphasizes the synchronization of the LFP segments around the spike times.

If all the spikes in bursts are used to calculate the coherence, the SFC values will decrease, even if there is strong phase-locked synchronization between bursty spike trains and LFPs at high frequency bands. To overcome this drawback, we propose an improvement to the SFC algorithm. The manipulation is analogous to the process of weighting, thus the modified

approach is referred to as weighted spike field correlation (WSFC). The weighting process emphasizes the contributions of some aspects of a dataset to the final result by allocating them greater weight during the analysis. When calculating the WSFC, only the first spike in each burst is employed to compute the SFC. The number of times that the first spike is used is determined by the number of spikes per burst. This manipulation emphasizes the firing time of the first spike in a burst and highlights the difference between a burst and an individual spike. To evaluate the performance of the proposed method, we applied it to simulation data and real neurobiological signals recorded in the hippocampus of rats.

2. Materials and Methods

2.1. Methods

The SFC is a function of frequency, which is obtained by computing the ratio of the power spectrum for the spike-triggered average (STA) over the average power spectrum of the LFP fractions (Fries, et al., 1997). Thus, the SFC is dependent on the LFP power and the spike number. Suppose that the spike train of a neuron is denoted as S   s1 ,s2 ,

sm  , where m is

the spike number. V  v1 ,v2 , vm  is the set of LFP segments, where vi is the sample of the LFP signal in the time window  si  T 2 , si  T 2 , and T is the duration of the LFP segments. The STA is constructed by averaging the LFP fractions within the windows centered on the spikes. The power spectrum of STA (PSTA) is defined as:

1

m



    vi  ,  m i 1 

(1)

where  denotes the operation used to calculate the power spectrum. To describe the power of every frequency component in the LFP segments used to construct the STA, i.e., vi with

i  1, 2,

,m , the average power spectrum of vi is:



1 m    vi  . m i 1 

(2)

This is also referred to as the spike-triggered power spectrum or STP (Fries, et al., 1997; Rutishauser, et al., 2010), and the SFC is defined as (Fries, et al., 1997):

 SFC 

 100% . 

(3)

The STP and PSTA can be computed using many methods. For example, multitaper analysis is a powerful and robust method for estimating a single-trial spectrum (Jarvis & Mitra, 2001), which can be performed using the Chronux toolbox (Bokil, Andrews, Kulkarni, Mehta, & Mitra, 2010). The multitaper method is employed in the present study to analyze the spectra of the simulated and experimentally recorded LFP signals.

The SFC reflects the synchronization between spike trains and LFPs at different frequencies. However, it does not function well with bursty spike trains and LFPs at high frequency bands, as shown in the following section. A burst can be defined as a temporary increase in the firing rate of spikes relative to the background activity (Cocatre-Zilgien & Delcomyn, 1992; Palm, 1981; Robin, et al., 2009). As the mechanism for generating bursts has been mentioned, it is commonly accepted that small depolarization keeps the cell silent, moderate depolarization makes the cell fire single spikes, and high depolarization causes the cell to discharge in burst mode (Harris, Hirase, Leinekugel, Henze, & Buzsáki, 2001). Thus, bursts code the same neural information as single spikes but with higher reliability (Harris, et al., 2001; Lisman, 1997). Based on this concept, the first spike in each burst is selected and used to represent the burst as an event (Kepecs & Lisman, 2003; Swadlow & Gusev, 2001).

The present study focuses on the approach used to quantify the level of synchronization between spike trains and LPFs. Thus, we propose the WSFC method, which allows only the first spike in each burst to be utilized in the computation of the PSTA and STP. To emphasize the difference between single spikes and bursts, the first spikes in bursts are used multiple times. The number of times is equal to the number of spikes per burst. In this manner, the first spike timing represents the occurrence of the burst and the weighting procedure (using the first spike multiple times) reflects the properties of the burst. Suppose that a bursty spike train is SB   s1 ,

, sm , b11 ,

, b1k1 , b21 ,

, b2k2 ,

, bn1 ,

, bnkn  , where s denotes single spikes, b

represents bursts, m is the number of single spikes, n is the number of bursts, and k is the

spike number per burst. The power spectrum of the weighted STA is 

n  m  1   vi   k j  v j1   , j 1   m   k  i 1

 w   

(4)

and the weighted spike triggered power spectrum of LFP is

w 

n  m  1    vi   k j  v j1  , m   k  i 1 j 1 

(5)

where v are the samples of the LFP signal around the spikes. Thus, the WSFC is defined as:

WSFC 

w 100% . w

(6)

In the present study, “weighted” means that the LFP segments around the first spikes in bursts are reused for calculating the STA and STP, where the “weight” is equal to the spike number per burst. The main difference between these two methods is that WSFC uses the first spike in a burst multiple times, whereas all of the spikes are used by SFC. This operation highlights the difference between bursts and single spikes. An important advantage is that our approach can remove the influence of second and later spikes in bursts during the computation of WSFC, which yields a more reliable estimate of the correlation between bursty spike trains and LFPs.

Similar to the SFC, the WSFC is a population method that cannot be calculated for single spikes. It is also affected by spike count used in the computation. Thus, it is necessary to minimize or avoid the bias generated by the amount of spikes, where two measures may be adopted. First, sufficient spikes (generally > 50) must be used to calculate the WSFC. Second, it is necessary to ensure that the spike counts are equal in different conditions (Rutishauser, et al., 2010). In addition, two parameters need to be identified before calculating the SFC and WSFC. The traces of the LFP used to construct the STA were set to 960 ms in the simulations, which facilitated the examination of low frequencies. For the real data used in the present study, a relatively short window of 480 ms was employed to focus on the gamma band frequencies. To summarize, the window length is selected according to the following principle: it should be sufficiently long to allow a reliable estimate of the power spectrum but short enough to represent the dynamics of the LFP signal in the desired band. In the simulations, the spectrum was estimated using the multitaper method with a time-bandwidth

product of four and seven tapers, which yielded a spectral resolution (half-bandwidth) of 4.2 Hz. During the analysis of real data, due to the short window of 480 ms, we used a time-bandwidth product of three and five tapers to obtain a relatively smaller spectral resolution of 6.25 Hz. This means that we need to select an appropriate time-bandwidth product and taper count to balance the spectral resolution and the benefit of the spectral estimate.

2.2. Simulated data

A superimposed waveform with two sinusoidal oscillations is used to demonstrate the calculation and interpretation of the SFC and WSFC (Fries, et al., 2002). Two sinusoids of 50-s duration are used for illustrative purposes, where one has a frequency of 10 Hz with an amplitude of 5 μV and the other has a frequency of 25 Hz with an amplitude of 2 μV , as shown in Fig. 1(a) and 1(b). The composite waveform is shown in Fig. 1(c). The rhombuses in these three plots indicate the occurrence of simulated spikes, which are phase-locked to a random subset of the troughs of the 25-Hz oscillation. At the same time, they are not synchronized with the 10-Hz component. This model with 500 spikes or bursts was used to generate the results presented in Section 3.1.

To simulate the real extracellular recoding, the LFP signal was generated by summing multiple sine waves with different frequencies, amplitudes, and phases (Rutishauser, et al., 2010). The frequencies ranged from 1 Hz to 100 Hz with a step of 1 Hz, thereby focusing on the LFP in the gamma band and below. The amplitudes of the components are inversely proportional to their frequencies. The phases were selected randomly from

0 2  .

In

addition, a white Gaussian noise with the signal-to-noise ratio of –10 dB was added to the composite oscillation. Thus, the artificial LFP signal generated followed the 1/f power distribution. The phase-locked spikes that are fired by simulated individual neurons are located at a certain phase of the underlying oscillation and they skip cycles at random, where the amount is denoted by n p , whereas the non-phase-locked spikes occur randomly and their

amount is denoted by nr . The ratio R  n p

n

p

 nr  determines the strength of the

correlation between spikes and LFPs, e.g., R  1 implies perfect phase locking and R  0 indicates complete non-phase locking. This model was used to generate the results presented in Section 3.2 and 3.3.

In these two models, the simulated bursts comprised 2–6 spikes. The inter-spike intervals (ISIs) in the bursts ranged from 3 to 10 ms. The amount of bursts was quantified using the burst index, which is defined as the ratio of spikes in bursts relative to all spikes (Mizuseki, Diba, Pastalkova, & Buzsaki, 2011; Mizuseki, Royer, Diba, & Buzsaki, 2012).

2.3. Experimentally recorded data

The dataset was recorded at the Gyorgy Buzsáki laboratory, New York University, and downloaded from the CRCNS Data Sharing website. This dataset has been used to analyze the relationship between theta oscillations and neuronal firing in the entorhinal-hippocampal system (Mizuseki, et al., 2009). The dataset comprises simultaneous recordings of cells in the CA1 layer of the right dorsal hippocampus of three Long-Evans rats, which were implanted with 4-shank or 8-shank silicon probes. After recovery from surgery (about one week) and training ( 3 days), the physiological signals were recorded during open field tasks, where the rats chased randomly dispersed drops of water or pieces of “Froot Loops” on an elevated square platform. The signals were amplified (1000), band-pass filtered (1–5000 Hz), and sampled at 20 kHz. The LFPs were then downsampled to 1250 Hz. Offline spike sorting was performed automatically. Further details of the behavioral experiment and data collection have been published previously (Mizuseki, et al., 2009). A subset of this dataset was analyzed in the present study (specifically, ec013.527 at http://crcns.org/data-sets/hc/hc-2).

3. Results

3.1. Demonstration of SFC and WSFC

In this section, we demonstrate the utilization of SFC and WSFC, and explain why the latter performs better with bursty spike trains. First, we illustrate the calculation of the SFC with perfect phase locking. The temporal waveform (±200 ms around the spike occurrence) and power spectrum of the STA are shown in Fig. 1(d) and Fig. 1(e), respectively. There is perfect synchronization between the spikes and the 25-Hz component, so this component was almost not attenuated by the averaging procedure. By contrast, the 10-Hz component is not visible because it was not consistently phase-locked to the spikes, thus it was averaged out. Figure 1(f) shows a plot of the average power spectrum of the LFP segments (STP). The SFC obtained is shown in Fig. 1(g). The value at a frequency of 25 Hz is 100%, thereby indicating that all of the spikes appear at exactly the same phase relative to the LFP component at this frequency. By contrast, the value at a frequency of 10 Hz is close to 0%, which means that the spikes lack phase synchronization with this frequency component. The non-zero values, except 100%, were caused mainly by insufficient averaging due to the finite spike number and the error of the fast Fourier transform (FFT) algorithm in MATLAB. In brief, the process used to calculate the STA and normalize the power spectrum of the STA determines that the SFC is dependent on the spike number and the power spectrum of the LFP. The SFC ranged from 0% to 100% and it varied as a function of the frequency. This value can describe the accuracy with which the spikes follow a particular phase of any frequency components in the LFP.

Next, some spikes were added to the simulated spike train and their effects on the SFC were investigated. These spikes were located at the peaks of the 25-Hz oscillation, which appear immediately after the spikes denoted by the rhombuses in Figs. 1(a–c), where they are indicated by the asterisks. In this case, the STA, PSTA, STP, and SFC are shown in Figs. 1(h–k). In these figures, the main frequency component is the 10-Hz oscillation. However, the amplitude of this component is severely reduced (about 0.02 μV ). The waveforms of the 25-Hz component in the windows around the spikes indicated by the rhombuses and the asterisks have the same amplitude but the opposite polarity, which makes the 25-Hz component disappear completely. However, there are two peaks at frequencies of 10 Hz and 25 Hz in the STP, which correspond to the two components in the LFP. Due to these two factors, the SFC values at all frequencies are approximately 0%. Therefore, we may conclude

that the spikes at a half-period of the frequency component after the phase-locked spikes (any phase and not necessarily in the trough) can degrade the SFC performance significantly. In the present study, these spikes are referred to as negative spikes because of their negative effect on the SFC value. The negative spikes could appear in bursts in the experimentally recorded spike trains, which is discussed in more detail later. In addition, it should be noted that the spikes that occurred exactly at the peaks were considered in this example, but the spikes around the peaks could become negative spikes and decrease the SFC value. We assumed that the negative spikes around the peaks were distributed uniformly within a window of 10-ms centered on specific peaks and the influence of negative spikes on the SFC was investigated. In Fig. 2, the simulation results show that these two types of negative spikes had similar effects on the calculation of the SFC. The SFC value declined by more than 50% when the number of negative spikes reached 25% of the number of phase-locked spikes in the simulated spike train. Furthermore, the SFC declined by more than 80% when the number of negative spikes comprised 100% of the phase-locked spikes.

8 4 0 -4 -8 8 4 0 -4 -8 8 4 0

5

5.1

2

0

2

5.6

5.7

5.8

5.9

6

(g) 100

10 5

50

(i)

0

0 Time (s)

0.2

0 0

10 20 30 40 50 Frequency (Hz) (j)

0.1

Power ( μV 2 )

0

-0.1 -0.2

0

0.2

0

0

10 20 30 40 50 Frequency (Hz)

0 0

10 20 30 40 50 Frequency (Hz) (k)

15

10 20 30 40 50 Frequency (Hz)

5

SFC (%)

0 Time (s)

0.1

Amplitude (mV)

5.4

-2 -0.2

(h)

(e) 4

5.5 Time (s) (f) 15

5.3

Power ( μV 2 )

Amplitude (mV)

(d)

5.2

SFC (%)

-8

Power ( μV 2 )

-4

Power ( μV 2 )

(c)

Amplitude (mV)

(b)

Amplitude (mV)

Amplitude (mV)

(a)

10

2.5

5 0 0

10 20 30 40 50 Frequency (Hz)

0

0

10 20 30 40 50 Frequency (Hz)

Fig. 1. Illustration of the calculation of the SFC. (a) Waveform of the 10-Hz component with an amplitude of 5 μV . (b) Waveform of the 25-Hz component with an amplitude of 2 μV . (c) Superimposed waveform of the two sinusoidal oscillations. The rhombuses in the three plots (a–c) indicate the occurrence of spikes that are phase-locked to a random subset of the troughs of the 25-Hz component. The asterisks denote the spikes that occur at peaks immediately after the rhombuses. (d) Waveform of the STA, which is located at ±200 ms around the spike occurrence. (e) Power spectrum of the STA. (f) Average power spectrum of the LFP segments used to construct the STA. (g) SFC between the simulated spikes and LFP. These four plots (d–g) were calculated using only the spikes marked by rhombuses. The other four plots (h–k) are the same as (d–g), but they were computed using all the simulated spikes.

1 a

100

0.81

around b exact a

75

SFC (%)

0.8 0.6 0.6 0.4

b

50 25

0.4 0.2

0 0

0.2 0 0

25

25 50

50 75 Percent75(%)

100 100

0

0 25 75 100 Fig. 2. Effects of two types of negative spikes50on the SFC. The x-axis label “Percent” denotes the ratio of

negative spikes relative to phase-locked spikes. The green line with triangles, which are labeled as “exact,” shows the effect of spikes that occurred exactly at the peaks. The blue line with circles, which are labeled as “around,” shows the effect of spikes that occurred around the peaks.

Finally, we demonstrate the power of WSFC to estimate the correlation between bursty spike trains and LFPs. In Fig. 3, the bursts are randomly phase-locked to the troughs of a 25-Hz sine wave. The vertical lines denote the occurrence of spikes and the triangles indicate the first spikes in the bursts. At 25 Hz, the WSFC and SFC were 100% and 4.3%, respectively. Thus, the second and later spikes in each burst reduce the SFC value severely,

Amplitude ( mV)

even with perfect synchronization. By contrast, the WSFC reflects the actual synchronization.

2 0 -2 1

1.2

1.4

1.6

1.8

2

Time (s)

Fig. 3. Phase locking bursts and 25-Hz component of the simulated LFP within a 1-s time window.

3.2. Rate dependency of the SFC and WSFC

The frequency-dependent measure of the linear spike-field correlation is known to depend on the overall neural activity, i.e., the spike rate or total number of spikes (Lepage, Kramer, & Eden, 2011). This rate dependency makes it difficult to compare the spike-field correlation across different experimental conditions. A procedure based on a generalized linear model has been proposed for separating the effects of the overall neural activity from the spike train and LFP oscillation coupling (Lepage, et al., 2013). The SFC and WSFC methods are also based

on spectral analysis. In order to assess the effects of neural activity on the SFC and WSFC, this section considers the rate dependency of these two methods. The SFC was examined in three cases: all of the spikes were locked to a certain phase of the 50-Hz component in the simulated LFP signal (100% phase-locked), half of the spikes were

phase-locked

(50%

phase-locked),

or

the

spikes

were

fired

randomly

(non-phase-locked). The duration of the simulated signals was 10 s. The spike rates ranged from 1 to 30 spikes/s with a step of 1 spikes/s. Fig. 4(a) shows the mean value of the SFC with different spike rates, which are based on 100 realizations for each rate. The WSFC was also examined in three cases: all the bursts were locked (100% phase-locked), half of the bursts were locked (50% phase-locked), or the bursts occurred randomly (non-phase-locked). The burst rate varied from 0.25 to 5 bursts/s with a step of 0.25 bursts/s. It was assumed that the bursts were only in the simulated spike trains. The mean WSFC values based on 100 realizations are plotted in Fig. 4(b). It shows that either the SFC or the WSFC was strongly dependent on the spike rate or burst rate, and they decreased as the rate increased. However, they tended to be stable when the rate was sufficiently large. We also note that the relative difference between the various phase-locked conditions remained almost the same, regardless of the rates of the spikes or bursts. This is of crucial importance when making comparisons between different experimental conditions. Furthermore, the results are consistent with those obtained based on spike numbers in a previous study (Rutishauser, et al., 2010), which showed that the SFC remained unchanged if sufficient spikes were used in the calculation. 0.3

Indeed, the rate and number are equivalent to some extent. 0.3 a100% phase locking 0.2 b c

(a)0.20.3 0.3

a b cRandom phase locking

a b50% phase locking c (b) 0.4 0.4

0.3 0.3

0.1

SFC

WSFC

0.2 0.10.2

0

0.1 0.1

0 0.25

0

0.25

0.5

0.75

00

0

0.5

0

0.75

0.25

1

0.2 0.2

0.5

0.1 0.1

0.75

1

1

5

10 10

15 20 Spike rate

25

30 30

00 0

1

2 3 Burst rate

4

55

Fig. 4. Illustration of the rate dependencies of the two methods. (a) Dependency of the SFC on the spike rate. (b) Dependency of the WSFC on the burst rate.

3.3. Effects of bursts in spike trains on the SFC and WSFC Using the simulated LFP signal, we investigated the effects of bursts in spike trains on the SFC and WSFC calculations, which demonstrated the advantage of the proposed method for estimating the correlation between bursty spike trains and LFPs. The SFC and WSFC were calculated for three cases: all the spikes and bursts were locked to a certain phase of specified components (50 Hz and 5 Hz for example) in the simulated LFP signal (100% phase-locked), half the spikes and bursts were phase-locked (50% phase-locked), or the spikes and bursts were fired randomly (non-phase-locked). The mean values of the SFC and WSFC at the two frequencies versus the burst indexes are plotted in Figs. 5(a–d). There were 100 realizations for each index. To avoid any bias due to the number of spikes, 200 spikes or bursts were generated for each realization. For the 50-Hz component, the SFC was reduced greatly when the bursts occurred in spike trains in both the 100% and 50% phase-locked cases, whereas it remained at the chance level with random phase-locking, as shown in Fig. 5(a). This means that the bursts did not lead to a spurious increase in the SFC with random spikes, but the SFC values can be reduced severely if the spikes are phase-locked to high frequency components. For the 5-Hz component, a decrease in the SFC also occurred with the 100% and 50% phase-locked cases, but with a smaller amplitude. However, the SFC increased slightly with random phase locking, as shown in Fig. 5(c). This indicates that the bursts have a relatively small effect on the SFC estimate when the spikes are phase-locked to low frequency components. All of the simulations were implemented using 200 spikes, thus the decrease in the SFC was caused by the location rather than the amount of spikes. By contrast, the WSFC method based on the first spikes effectively eliminated the effects caused by bursts at both low and high frequencies. The reuse of the first spikes in bursts means that the amplitude of the STA increases as a function of the number of bursts involved. Thus, the WSFC exhibited only a slight increase as the bursts in spike trains increased, as shown in Fig. 5(b) and 5(d). However, the relative difference in the WSFC between conditions remained almost the same, which is useful when making comparisons between different phase-locked cases.

0.3 a b cRandom

0.3 a b50% c

100% phase locking a 0.2 b c

0.2

(a) 20

10

0 0.25

5

0.75

0.5

0

0.25 0.75

15 10

0.5 1

0.75

0 0

0.25

0.5 0.75 Burst Index

0

1

(c) 20

(d) 30

15

20

10 5 0

1

5

1

WSFC (%)

0.5

0

SFC (%)

0.25

0.1

WSFC (%)

SFC (%)

0

phase locking

(b) 20

15

0.1

phase locking

0

0.25

0.5 0.75 Burst Index

1

0

0.25

0.5 0.75 Burst Index

1

0

0.25

0.5 0.75 Burst Index

1

10 5 0

Fig. 5. (a), (b) SFC and WSFC between bursty spike trains and the LFP in different phase-locked cases. Spikes or bursts are phased locked to the 50-Hz component. (c), (d) SFC and WSFC between bursty spike trains and the LFP, where phase locking is shown between spikes or bursts and the 5-Hz component.

Furthermore, we examined the influence of bursts at the following nine frequencies: 5 Hz in the theta band; 10 Hz in the alpha band; 20 Hz in the beta band; and 30 Hz, 40 Hz, 50 Hz, 60 Hz, 70 Hz, and 80 Hz in the gamma band. Single spikes and the first spikes in bursts were assumed to fire exactly at a certain phase of these component oscillations. The effects were quantified based on the coefficient of variation (CV) of the SFC or WSFC magnitudes with different burst indexes. As shown in Fig. 6(a) and 6(b), the SFC between spike trains and LFPs in the gamma frequency band tended to be affected by bursts in the 100% and 50% phase-locked cases. However, with low frequencies, e.g., 5 Hz and 10 Hz, the bursts in the spike trains had very little effect on the computation of the SFC. By contrast, the WSFC exhibited relatively robust performance with both high and low frequencies. The effects of bursts with random phase-locking are shown in Fig. 6(c). The SFC and WSFC measures exhibited similar performance. The variation in the WSFC was relatively large, but this did not indicate a reduction in the WSFC performance. As explained earlier, these variations are present in different locked conditions with similar increasing trends, thereby allowing reliable comparisons of the correlations between different conditions. This was also demonstrated by

the CV values in the different phase-locked cases. With a lower degree of phase locking, the mean of the WSFC was smaller at different burst indexes and the CV was higher at the frequencies considered. In summary, the WSFC is an efficient tool for determining the correlations between spike trains and LFPs, particularly in high frequency bands.

0.6

0.6

0.4

30.4

CV

(b) 0.8

CV

(a) 0.8

0.2

0.2 0 (c)

2.5 5

10 20 30 40 50 60 70 80 Frequency (Hz)

0.8

CV

0.6

0

2 1.5

0.4

2.5 5

10 20 30 40 50 60 70 80 Frequency (Hz)

2

1.5 SFC

0.2 0

3

5

10 20 30 40 50 60 70 80 Frequency (Hz)

1

1

0.5

0.5

WSFC

Fig. 6. Effects of bursts on the SFC and WSFC, where the spikes or bursts are phase-locked to different

0

0

0 2locking. 4 0 (c) 6 2 Random 8 4 6phase 8 locking. frequency components. (a) 100% phase locking. (b) 50% phase

3.4. Application to real data

In the hippocampus, gamma frequency oscillations (30–80 Hz) have been suggested to underlie various cognitive functions, such as attention selection (Bauer, Oostenveld, Peeters, & Fries, 2006; Fries, et al., 2001), memory (Fell, et al., 2001; Howard, et al., 2003; van Vugt, Schulze-Bonhage, Litt, Brandt, & Kahana, 2010), and sensory perception (Colgin, et al., 2009; Muzzio, et al., 2009). It has also been reported that the firing patterns of pyramidal cells in the hippocampus are significantly phase-locked to gamma oscillations in rats that exhibit behaviors (Colgin, et al., 2009; Csicsvari, et al., 2003; Senior, Huxter, Allen, O'Neill, & Csicsvari, 2008). Thus, we examined whether this phase locking can be characterized using the WSFC and SFC methods. The pyramidal cells exhibited firing patterns with single spikes and complex spike bursts in the recordings that comprised the dataset used in the present study. A segment (0–450 s) of the recordings used for correlation analysis is shown in Fig.

7(a) and 7(b), which contained the gamma band LFP signal and neuronal activity (neuron 37 in the selected dataset).

To preserve the timing relationship between spikes and the LFP, gamma bandpass filtering was performed digitally with a zero-phase shift using the EEGLAB toolbox (Delorme & Makeig, 2004). We then estimated the WSFC and SFC as a function of time (sliding window of 10 s advanced in steps of 5 s). A series of two or more consecutive spikes with ISI 50 spikes. (d) Z-transformed WSFC for each selected window. The red horizontal line indicates the significant at p  0.05 . (e) Z-transformed SFC for each selected window. The correlations in plots (d) and (e) were obtained using the average of the WSFC and SFC in the gamma band, respectively.

4. Discussion

LFPs comprise low frequency extracellular voltage fluctuations, which reflect the neural mass activity. Spike trains contain series of action potentials that are fired by individual neurons and they can describe the activity of a single neuron. Determining the relationship between spiking activities and LFPs is an important tool for investigating the rhythmic neuronal synchronization of brain functions. The SFC is used extensively as a measure to

estimate the correlation between spike trains and LFPs. Applications of the SFC method to simulated and experimental data have demonstrated that the bursts in spike trains degrade its performance in proportion to the burst index. Some neurons fire without bursts or with small bursts, but it cannot be assumed that bursts never occur. Indeed, various neuronal types exhibit typical bursting behaviors in the central nervous system (Connors, Gutnick, & Prince, 1982; Goldberg, et al., 2002; Gray & McCormick, 1996; Harris, et al., 2001). In these cases, the bursts in spike trains have significant effects on the estimated SFC. To address this issue, we propose a correction method called WSFC, which utilizes the first spikes multiple times, rather than the original spikes in bursts. The number of times that the first spike is used is equal to the spike count per burst. The concept that underlies this manipulation is as follows. Bursts are produced either by strong synaptic stimulations or by depolarizing current injections. As a first approximation, the instantaneous firing rate of neurons is proportional to the magnitude of the input (Kepecs & Lisman, 2003). Previous studies have shown that the timing of the first burst spike is correlated with the stimulus features (Lesica & Stanley, 2004; Oswald, Chacron, Doiron, Bastian, & Maler, 2004). Based on these considerations, multiple utilizations of the first spike are employed to preserve the occurrence and duration information for each burst. The main aim of this manipulation is to remove the effects of the second and latter spikes in bursts on the computation of the SFC. This method also aims to ensure that the number of spikes remains unchanged, thereby eliminating the bias generated by different amounts of spikes.

The proposed correction method was shown to produce major improvements in the correlation estimate, where the WSFC did not decrease with increases in the burst index; indeed it actually increased slightly. The WSFC varied as a function of the burst index, but the relative difference in the WSFC between different phase-locked conditions remained almost the same, thereby allowing reliable comparisons of the degree of phase synchronization. Thus, the WSFC is more suitable than the SFC for characterizing the phase locking between bursty spike trains and LFPs. Next, we explain the mechanism that underlies the increase in the WSFC. The reuse of the first spike in a burst allows the same segment of the LFP to be employed in the computation of the STA many times, so the amplitude of the STA increases, where the increase is proportional to the burst index. As a result, the PSTA

increases with the burst index and the WSFC increases as a consequence.

We also note that the SFC increases and it does not decrease when the spikes are randomly phase-locked to the LFP component at low frequencies (e.g., 5 Hz in Fig. 5(b)). The reason for this contrasting result is explained below. In the low frequency band, the period of the LFP is much longer than that of the ISIs in bursts, thus the segments of the LFP centered on each spike in a burst are similar to each other. Assuming that there is no phase locking, this makes the amplitude of the STA increase and the SFC becomes larger as a function of the bursts considered. By contrast, due to the short period of the high frequency LFP component and random phase locking, the amplitude of the STA decreases and the SFC becomes smaller when sufficient spikes are involved.

Previous studies have shown that pyramidal cells in the hippocampus are phase-locked to gamma oscillations in rats that exhibit behaviors (Colgin, et al., 2009; Csicsvari, et al., 2003; Senior, et al., 2008). This phase relationship is described by the phase distribution of the spikes fired by a neuron. The phase of each spike is assigned as the phase of the normalized gamma cycle at the spiking time. These gamma epochs are determined based on the mean and standard deviation of the recording signal. However, this detection method is not capable of detecting less regular or lower amplitude gamma oscillations (Senior, et al., 2008), which may result in a biased estimate of the phase distribution and, consequently, the coupling between the unit activity and gamma oscillations. However, the WSFC method is not affected by this issue and it can effectively characterize the degree of phase locking at different frequencies. The WSFC behaves in a similar manner to the SFC when the spike phase is locked to the low frequency components of LFP, but it performs much better with bursty spike trains and high frequency components.

In the present study, we defined a burst as a short period with high firing rates (Mammone & Morabito, 2008). However, different temporal structures exist in bursts in various species, which play different roles in information encoding (Kepecs & Lisman, 2003; Metzner, Koch, Wessel, & Gabbiani, 1998). In the electric fish, for example, the ISIs within bursts might be as long as 25 ms (Metzner, et al., 1998), whereas the typical ISIs range from 4 to 10 ms in

cortical and hippocampal neurons (Harris, et al., 2001; Reich, Mechler, Purpura, & Victor, 2000). Bursts have different temporal structures, but they are triggered when the input current of the neurons increases and the spike count per burst encodes the degree of the increase (Kepecs, Wang, & Lisman, 2002). When calculating the WSFC, the weight can reflect this encoding process because it is equal to the spike count per burst. Thus, it can be inferred that the WSFC will work well for all types of bursts. Before calculating the WSFC, it is important to detect the bursts in spike trains using either a fixed threshold or self-adapting methods (Gourevitch & Eggermont, 2007; Pasquale, Martinoia, & Chiappalone, 2010; Tokdar, Xi, Kelly, & Kass, 2010). During artifact rejection from electroencephalographic recordings, the procedure based on the estimation of kurtosis and Renyi‟s entropy can detect the artifacts automatically (Mammone & Morabito, 2008). This method could be used to detect bursts in spike trains and our future research will address the improvement of automatic burst detection.

5. Conclusion

Bursts are often present in recorded spike trains in real conditions and they are believed to contain more reliable information than single spikes. However, the widely used SFC method underestimates the correlation between bursty spike trains and LFPs, especially in the high frequency bands. In the present study, we proposed a correction measure called WSFC, which allows only the first spike in each burst to be utilized in the calculation procedure, thereby eliminating the effects of the other spikes. WSFC also utilizes the first spike multiple times to highlight the difference between single spikes and bursts. This „weighting‟ operation preserves the inherent characteristic of bursts, i.e., a series of spikes generated by a large depolarization. The WSFC is a rate-dependent measure, but it is capable of allowing reliable comparisons between different experimental conditions by using large numbers of spikes (bursts) or by taking the same number of spikes (bursts) across conditions. Our simulation and experimental results showed that the proposed method performed better than SFC when investigating the relationship between bursty spike trains and high frequency band LFPs. It can also be used to analyze any type of spike train and LFP in any frequency band. Furthermore, the WSFC could potentially be applied to study whether bursts enhance the

correlation between spike trains and LFPs. In summary, the WSFC is a promising method for elucidating the details of neural coding.

Acknowledgments This research was supported in part by the National Natural Science Foundation of China (61025019, 61105027 and 61203210), and the Research and Development Program of Science and Technology in Qinhuangdao, China (2012021A047). We thank the Gyorgy Buzsáki laboratory for collecting the data and making it available through the CRCNS program (http://crcns.org).

References Bauer, M., Oostenveld, R., Peeters, M., & Fries, P. (2006). Tactile spatial attention enhances gamma-band activity in somatosensory cortex and reduces low-frequency activity in parieto-occipital areas. Journal of Neuroscience, 26, 490-501. Bokil, H., Andrews, P., Kulkarni, J. E., Mehta, S., & Mitra, P. P. (2010). Chronux: a platform for analyzing neural signals. Journal of Neuroscience Methods, 192, 146-151. Chalk, M., Herrero, J. L., Gieselmann, M. A., Delicato, L. S., Gotthardt, S., & Thiele, A. (2010). Attention reduces stimulus-driven gamma frequency oscillations and spike field coherence in V1. Neuron, 66, 114-125. Claverol-Tinture, E., Cabestany, J., & Rosell, X. (2007). Multisite recording of extracellular potentials produced by microchannel-confined neurons in-vitro. IEEE Transactions on Biomedical Engineering, 54, 331-335. Cocatre-Zilgien, J. H., & Delcomyn, F. (1992). Identification of bursts in spike trains. Journal of Neuroscience Methods, 41, 19-30. Colgin, L. L., Denninger, T., Fyhn, M., Hafting, T., Bonnevie, T., Jensen, O., Moser, M. B., & Moser, E. I. (2009). Frequency of gamma oscillations routes flow of information in the hippocampus. Nature, 462, 353-357. Connors, B. W., Gutnick, M. J., & Prince, D. A. (1982). Electrophysiological properties of neocortical neurons in vitro. Journal of Neurophysiology, 48, 1302-1320. Courtemanche, R., Pellerin, J. P., & Lamarre, Y. (2002). Local field potential oscillations in primate cerebellar cortex: modulation during active and passive expectancy. Journal of Neurophysiology, 88, 771-782. Csicsvari, J., Jamieson, B., Wise, K. D., & Buzsaki, G. (2003). Mechanisms of gamma oscillations in the hippocampus of the behaving rat. Neuron, 37, 311-322. Delorme, A., & Makeig, S. (2004). EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis. Journal of Neuroscience Methods, 134, 9-21. Eggermont, J. J., & Smith, G. M. (1995). Synchrony between single-unit activity and local field potentials in relation to periodicity coding in primary auditory cortex. Journal of Neurophysiology, 73, 227-245. Fell, J., Klaver, P., Lehnertz, K., Grunwald, T., Schaller, C., Elger, C. E., & Fernandez, G. (2001). Human memory formation is accompanied by rhinal-hippocampal coupling and decoupling. Nature Neuroscience, 4, 1259-1264. Fries, P., Reynolds, J. H., Rorie, A. E., & Desimone, R. (2001). Modulation of oscillatory neuronal synchronization by selective visual attention. Science, 291, 1560-1563. Fries, P., Roelfsema, P. R., Engel, A. K., Konig, P., & Singer, W. (1997). Synchronization of oscillatory responses in visual cortex correlates with perception in interocular rivalry. Proceedings of the National Academy of Sciences of the United States of America, 94, 12699-12704. Fries, P., Schroder, J. H., Roelfsema, P. R., Singer, W., & Engel, A. K. (2002). Oscillatory neuronal synchronization in primary visual cortex as a correlate of stimulus selection. Journal of Neuroscience, 22, 3739-3754. Galashan, F. O., Rempel, H. C., Meyer, A., Gruber-Dujardin, E., Kreiter, A. K., & Wegener, D. (2011). A new type of recording chamber with an easy-to-exchange microdrive array for chronic recordings in macaque monkeys. Journal of Neurophysiology, 105, 3092-3105. Goldberg, J. A., Boraud, T., Maraton, S., Haber, S. N., Vaadia, E., & Bergman, H. (2002). Enhanced synchrony among primary motor cortex neurons in the 1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine primate model of Parkinson's disease. Journal of Neuroscience, 22, 4639-4653. Gourevitch, B., & Eggermont, J. J. (2007). A nonparametric approach for detection of bursts in spike trains.

Journal of Neuroscience Methods, 160, 349-358. Grasse, D. W., & Moxon, K. A. (2010). Correcting the bias of spike field coherence estimators due to a finite number of spikes. Journal of Neurophysiology, 104, 548-558. Gray, C. M., & McCormick, D. A. (1996). Chattering cells: superficial pyramidal neurons contributing to the generation of synchronous oscillations in the visual cortex. Science, 274, 109-113. Gregoriou, G. G., Gotts, S. J., Zhou, H., & Desimone, R. (2009). High-Frequency, Long-Range Coupling Between Prefrontal and Visual Cortex During Attention. Science, 324, 1207-1210. Hagan, M. A., Dean, H. L., & Pesaran, B. (2012). Spike-field activity in parietal area LIP during coordinated reach and saccade movements. Journal of Neurophysiology, 107, 1275-1290. Harris, K. D., Henze, D. A., Hirase, H., Leinekugel, X., Dragoi, G., Czurko, A., & Buzsaki, G. (2002). Spike train dynamics predicts theta-related phase precession in hippocampal pyramidal cells. Nature, 417, 738-741. Harris, K. D., Hirase, H., Leinekugel, X., Henze, D. A., & Buzsáki, G. (2001). Temporal Interaction between Single Spikes and Complex Spike Bursts in Hippocampal Pyramidal Cells. Neuron, 32, 141-149. Howard, M. W., Rizzuto, D. S., Caplan, J. B., Madsen, J. R., Lisman, J., Aschenbrenner-Scheibe, R., Schulze-Bonhage, A., & Kahana, M. J. (2003). Gamma oscillations correlate with working memory load in humans. Cerebral Cortex, 13, 1369-1374. Issa, E. B., & Wang, X. (2011). Altered neural responses to sounds in primate primary auditory cortex during slow-wave sleep. Journal of Neuroscience, 31, 2965-2973. Jarvis, M. R., & Mitra, P. P. (2001). Sampling properties of the spectrum and coherency of sequences of action potentials. Neural Computation, 13, 717-749. Kepecs, A., & Lisman, J. (2003). Information encoding and computation with spikes and bursts. Network, 14, 103-118. Kepecs, A., Wang, X. J., & Lisman, J. (2002). Bursting neurons signal input slope. Journal of Neuroscience, 22, 9053-9062. Le Van Quyen, M., Bragin, A., Staba, R., Crepon, B., Wilson, C. L., & Engel, J., Jr. (2008). Cell type-specific firing during ripple oscillations in the hippocampal formation of humans. Journal of Neuroscience, 28, 6104-6110. Lee, H., Simpson, G. V., Logothetis, N. K., & Rainer, G. (2005). Phase locking of single neuron activity to theta oscillations during working memory in monkey extrastriate visual cortex. Neuron, 45, 147-156. Lepage, K. Q., Gregoriou, G. G., Kramer, M. A., Aoi, M., Gotts, S. J., Eden, U. T., & Desimone, R. (2013). A procedure for testing across-condition rhythmic spike-field association change. Journal of Neuroscience Methods, 213, 43-62. Lepage, K. Q., Kramer, M. A., & Eden, U. T. (2011). The dependence of spike field coherence on expected intensity. Neural Computation, 23, 2209-2241. Lesica, N. A., & Stanley, G. B. (2004). Encoding of natural scene movies by tonic and burst spikes in the lateral geniculate nucleus. Journal of Neuroscience, 24, 10731-10740. Lewandowski, B. C., & Schmidt, M. (2011). Short bouts of vocalization induce long-lasting fast gamma oscillations in a sensorimotor nucleus. Journal of Neuroscience, 31, 13936-13948. Lisman, J. E. (1997). Bursts as a unit of neural information: making unreliable synapses reliable. Trends in Neurosciences, 20, 38-43. Mammone, N., & Morabito, F. C. (2008). Enhanced automatic artifact detection based on independent component analysis and Renyi's entropy. Neural Networks, 21, 1029-1040. Metzner, W., Koch, C., Wessel, R., & Gabbiani, F. (1998). Feature extraction by burst-like spike patterns in multiple sensory maps. Journal of Neuroscience, 18, 2283-2300. Mizuseki, K., Diba, K., Pastalkova, E., & Buzsaki, G. (2011). Hippocampal CA1 pyramidal cells form

functionally distinct sublayers. Nature Neuroscience, 14, 1174-1181. Mizuseki, K., Royer, S., Diba, K., & Buzsaki, G. (2012). Activity dynamics and behavioral correlates of CA3 and CA1 hippocampal pyramidal neurons. Hippocampus, 22, 1659-1680. Mizuseki, K., Sirota, A., Pastalkova, E., & Buzsaki, G. (2009). Theta oscillations provide temporal windows for local circuit computation in the entorhinal-hippocampal loop. Neuron, 64, 267-280. Muzzio, I. A., Levita, L., Kulkarni, J., Monaco, J., Kentros, C., Stead, M., Abbott, L. F., & Kandel, E. R. (2009). Attention enhances the retrieval and stability of visuospatial and olfactory representations in the dorsal hippocampus. PLoS Biology, 7, e1000140. Oswald, A. M., Chacron, M. J., Doiron, B., Bastian, J., & Maler, L. (2004). Parallel processing of sensory input by bursts and isolated spikes. Journal of Neuroscience, 24, 4351-4362. Palm, G. (1981). Evidence, information, and surprise. Biological Cybernetics, 42, 57-68. Pasquale, V., Martinoia, S., & Chiappalone, M. (2010). A self-adapting approach for the detection of bursts and network bursts in neuronal cultures. Journal of Computational Neuroscience, 29, 213-229. Perelman, Y., & Ginosar, R. (2007). An integrated system for multichannel neuronal recording with spike/LFP separation, integrated A/D conversion and threshold detection. IEEE Transactions on Biomedical Engineering, 54, 130-137. Pesaran, B., Pezaris, J., Sahani, M., Mitra, P., & Andersen, R. (2002). Temporal structure in neuronal activity during working memory in macaque parietal cortex. Nature Neuroscience, 5, 805-811. Pienkowski, M., & Eggermont, J. J. (2011). Sound frequency representation in primary auditory cortex is level tolerant for moderately loud, complex sounds. Journal of Neurophysiology, 106, 1016-1027. Reich, D. S., Mechler, F., Purpura, K. P., & Victor, J. D. (2000). Interspike intervals, receptive fields, and information encoding in primary visual cortex. Journal of Neuroscience, 20, 1964-1974. Robin, K., Maurice, N., Degos, B., Deniau, J. M., Martinerie, J., & Pezard, L. (2009). Assessment of bursting activity and interspike intervals variability: a case study for methodological comparison. Journal of Neuroscience Methods, 179, 142-149. Rutishauser, U., Ross, I. B., Mamelak, A. N., & Schuman, E. M. (2010). Human memory strength is predicted by theta-frequency phase-locking of single neurons. Nature, 464, 903-907. Schwartz, A. B. (2004). Cortical neural prosthetics. Annual Review of Neuroscience, 27, 487-507. Senior, T. J., Huxter, J. R., Allen, K., O'Neill, J., & Csicsvari, J. (2008). Gamma oscillatory firing reveals distinct populations of pyramidal cells in the CA1 region of the hippocampus. Journal of Neuroscience, 28, 2274-2286. Siapas, A. G., Lubenov, E. V., & Wilson, M. A. (2005). Prefrontal phase locking to hippocampal theta oscillations. Neuron, 46, 141-151. Siegel, M., & Konig, P. (2003). A functional gamma-band defined by stimulus-dependent synchronization in area 18 of awake behaving cats. Journal of Neuroscience, 23, 4251-4260. Sirota, A., Montgomery, S., Fujisawa, S., Isomura, Y., Zugaro, M., & Buzsaki, G. (2008). Entrainment of neocortical neurons and gamma oscillations by the hippocampal theta rhythm. Neuron, 60, 683-697. Stafford, B. K., Sher, A., Litke, A. M., & Feldheim, D. A. (2009). Spatial-temporal patterns of retinal waves underlying activity-dependent refinement of retinofugal projections. Neuron, 64, 200-212. Swadlow, H. A., & Gusev, A. G. (2001). The impact of 'bursting' thalamic impulses at a neocortical synapse. Nature Neuroscience, 4, 402-408. Tiesinga, P. H., Fellous, J. M., Salinas, E., Jose, J. V., & Sejnowski, T. J. (2004). Inhibitory synchrony as a mechanism for attentional gain modulation. J Physiology Paris, 98, 296-314. Tokdar, S., Xi, P., Kelly, R. C., & Kass, R. E. (2010). Detection of bursts in extracellular spike trains using hidden semi-Markov point process models. Journal of Computational Neuroscience, 29, 203-212. van Vugt, M. K., Schulze-Bonhage, A., Litt, B., Brandt, A., & Kahana, M. J. (2010). Hippocampal gamma

oscillations increase with memory load. Journal of Neuroscience, 30, 2694-2699. van Wingerden, M., Vinck, M., Lankelma, J. V., & Pennartz, C. M. (2010). Learning-associated gamma-band phase-locking of action-outcome selective neurons in orbitofrontal cortex. Journal of Neuroscience, 30, 10025-10038. Vinck, M., Battaglia, F. P., Womelsdorf, T., & Pennartz, C. (2012). Improved measures of phase-coupling between spikes and the Local Field Potential. Journal of Computational Neuroscience, 33, 53-75. Vinck, M., Lima, B., Womelsdorf, T., Oostenveld, R., Singer, W., Neuenschwander, S., & Fries, P. (2010). Gamma-phase shifting in awake monkey visual cortex. Journal of Neuroscience, 30, 1250-1257. Vinck, M., van Wingerden, M., Womelsdorf, T., Fries, P., & Pennartz, C. M. (2010). The pairwise phase consistency: a bias-free measure of rhythmic neuronal synchronization. Neuroimage, 51, 112-122. Wang, Y., Iliescu, B. F., Ma, J., Josic, K., & Dragoi, V. (2011). Adaptive changes in neuronal synchronization in macaque V4. Journal of Neuroscience, 31, 13204-13213. Xu, S., Jiang, W., Poo, M. M., & Dan, Y. (2012). Activity recall in a visual cortical ensemble. Nature Neuroscience, 15, 449-455.