Subthalamic nucleus long-range synchronization—an independent hallmark of human Parkinson's disease

Beta-band synchronous oscillations in the dorsolateral region of the subthalamic nucleus (STN) of human patients with Parkinson's disease (PD) have been frequently reported. However, the correlation between STN oscillations and synchronization has not been thoroughly explored. The simultaneous recordings of 2390 multi-unit pairs recorded by two parallel microelectrodes (separated by fixed distance of 2 mm, n = 72 trajectories with two electrode tracks >4 mm STN span) in 57 PD patients undergoing STN deep brain stimulation surgery were analyzed. Automatic procedures were utilized to divide the STN into dorsolateral oscillatory and ventromedial non-oscillatory regions, and to quantify the intensity of STN oscillations and synchronicity. Finally, the synchronicity of simultaneously vs. non-simultaneously recorded pairs were compared using a shuffling procedure. Synchronization was observed predominately in the beta range and only between multi-unit pairs in the dorsolateral oscillatory region (n = 615). In paired recordings between sites in the dorsolateral and ventromedial (n = 548) and ventromedial-ventromedial region pairs (n = 1227), no synchronization was observed. Oscillation and synchronicity intensity decline along the STN dorsolateral-ventromedial axis suggesting a fuzzy border between the STN regions. Synchronization strength was significantly correlated to the oscillation power, but synchronization was no longer observed following shuffling. We conclude that STN long-range beta oscillatory synchronization is due to increased neuronal coupling in the Parkinsonian brain and does not merely reflect the outcome of oscillations at similar frequency. The neural synchronization in the dorsolateral (probably the motor domain) STN probably augments the pathological changes in firing rate and patterns of subthalamic neurons in PD patients.

Conclusive evidence of the correlation (and causality) between neuronal oscillations and synchronization in the PD STN has remained elusive. Physiological studies of neuronal synchronization in the STN of the MPTP primate model are not yet reported. Robust oscillatory synchronization patterns of STN spiking activity have been reported in the 6-hydroxydopamine rodent model of Parkinsonism (Machado et al., 2006;Mallet et al., 2008a,b;Lintas et al., 2012). In human PD patients, oscillatory synchronization of spiking activity has been reported in several studies (Levy et al., 2000(Levy et al., , 2002aAmirnovin et al., 2004;Weinberger et al., 2006;Hanson et al., 2012;Alavi et al., 2013;Lourens et al., 2013) but there have been no detailed descriptions of the dependence of the neuronal synchronization on the oscillatory activity or the spatial properties of the neuronal pairs (e.g., simultaneous recording of neurons from the oscillatory and non-oscillatory regions of the STN, see below).
Previous studies have shown that the STN of PD patients can be divided into a dorso-lateral oscillatory region (DLOR) and ventro-medial non-oscillatory region (VMNR) (Moran et al., 2008;Zaidel et al., 2010;Seifried et al., 2012;Guo et al., 2013). The first aim of this study was to explore the properties of neuronal (spike) synchronization of the STN of human PD patients, principally within and between the different STN domains. The second goal was to further explore the relationship between oscillations and synchronization phenomena in the neural activity of the STN.
To overcome the inherent technical difficulties of spike isolation (Joshua et al., 2007;Hill et al., 2011) and spike sorting (Lewicki, 1998) in the electrically noisy environment of the human operating room, and to increase the sensitivity of correlation analysis (Bedenbaugh and Gerstein, 1997;Gerstein, 2000) this study used the unresolved collective (multi-unit) spiking activity recorded by two different microelectrodes exploring the boundaries and the domains of the STN during DBS procedures. This enabled the exploration of the properties of long-range correlation in the STN, in contrast to correlation studies of the activity recorded by a single electrode (e.g., Moran et al., 2008) which can only reveal short range correlations.

PATIENTS AND SURGERY
Simultaneous microelectrode recordings from two electrodes in patients with Parkinson's disease (PD) undergoing surgery for subthalamic nucleus (STN) deep brain stimulation (DBS) were analyzed in this study. All patients met accepted criteria for STN DBS and signed informed consent for surgery. Microelectrode recording is performed to accurately localize STN borders and domains, in order to optimize the placement of the DBS electrode and thus enhance the therapeutic effects of the DBS procedure. The data collection was therefore done as part of our routine procedures and not part of a clinical trial. This study was authorized and approved by the Institutional Review Board of Hadassah University Hospital in accordance with the Helsinki Declaration (reference codes: 0545-08-HMO and HMO: 10-18.01.08).
Surgery was performed using a CRW stereotactic frame (Radionics, Burlington, MA, USA). STN target coordinates were chosen as a composite of the indirect anterior commissureposterior commissure (AC-PC) atlas-based location and direct (1.5 or 3 Tesla) T2 magnetic resonance imaging (MRI), using Framelink 4 or 5 software (Medtronic, Minneapolis, USA). The recordings used in this study were made while the patients were awake without sedation. The patient's level of awareness was continuously assessed clinically and, if drowsy, the patient was stimulated and awoken through conversation by a member of the surgical team. Data were obtained while the patients were off dopaminergic medication, which was stopped 12 h prior to surgery.

MICROELECTRODE RECORDINGS
Data were acquired with the MicroGuide system (Alpha-Omega Engineering, Nazareth, Israel). Neurophysiological activity was recorded using polyamide coated tungsten microelectrodes (Alpha Omega) with impedance mean ± standard deviation (SD) of 0.60 ± 0.11 M (measured at 1 kHz at the beginning of each trajectory). The signal was amplified by 10,000, band-passed filtered from 250 to 6000 Hz using four-pole Butterworth filter hardware, and sampled at 48 kHz by a 12-bit A/D converter (using ±5 V input range). Local field potentials were not recorded due to constraints of electrical noise in the operating room.
Microelectrode recording was performed using two parallel microelectrodes starting 10 mm above the estimated center of the dorsolateral STN target, based on the pre-operative T2 MRI image. The two electrodes were simultaneously advanced, and therefore the distance between the two electrodes was fixed (2 mm) during all recordings. Trajectories followed a doubleoblique approach (approximately 60 • from the axial AC-PC plane and 15 • from the mid-sagittal plane) toward the STN target. The angles of the trajectory were slightly modified to avoid the cortical sulci, the ventricles and major blood vessels as revealed by gadolinium-enhanced T1 MRI (Machado et al., 2006). The "central" electrode was directed at the center of the STN target, and an "anterior" (ventral) electrode was located 2 mm anterior to the central electrode. Typically, the electrodes were advanced in steps of ∼100μm between successive recordings sites within the STN. Only trajectories where both electrodes had passed through the STN for at least 4 mm were used in this study (yielding 72 trajectories of 2 electrodes from 57 PD patients undergoing bilateral STN deep brain stimulation surgery). After identification of the STN ventral border by the electro-physiologist, the STN and its sub-regions were automatically detected using the Hidden Markov model (HMM) method (Zaidel et al., 2009).

DATABASE
We studied 72 STN trajectories (each of 2 electrodes) from 57 PD patients, 40 males and 17 females, aged 58.9 ± 10.3 years (mean ± standard deviation, SD) and with disease duration of 10.3 ± 4.7 years (mean ± SD). The UPDRS motor part score, UPDRS III, was 49.2 ± 17.8 (mean ± SD) when assessed off dopamine replacement therapy before surgery. Patient details and clinical effects of the surgery are given at Table 1.
The minimal recording time duration of a STN pair to be included in this study was 5 s (analysis of the subset of recording with minimal recording duration of 10 s reveal similar results, data not shown). A total of 2390 multi-unit pair sites, in which both electrodes were judged to be inside the STN for the minimal   duration, were studied. The same data base was used for the single site (oscillation) analysis, yielding 4780 single STN sites. Recording (and analysis) time duration of the STN pairs equaled 23.7 ± 25.3 s (mean ± SD).

ANALYSIS OF SYNCHRONIZATION AND OSCILLATIONS
All data analysis utilized custom-made MATLAB 7.10b (R2010.b) routines. The local field potential frequency domain was filtered out by the recording apparatus. Burst frequencies below the range of the operating room band-pass filter (250-6000 Hz) could be detected using the rectified signal, which follows the envelope of multi-unit activity (Moran et al., 2008;Halliday and Farmer, 2010;Moran and Bar-Gad, 2010;Zaidel et al., 2010). The raw 250-6000 Hz analog signal was therefore rectified by the "absolute" operator and the global mean was subtracted. Thus, the resulting analysis represents only spike activity. The average power spectrum density (PSD) at each site was calculated using Welch's method with a 1.5 s Hamming window (50% overlap), after removing the local window mean, and with a 131,072 FFT points (nfft), yielding spectral resolution of 1/3 Hz [nfft = 2 ∧ round(log2(Fs/f_res)), where Fs = sampling frequency and f_res is the spectral resolution]. PSD amplitude is affected by the amplitude of the recorded neural activity, which is impacted by non-physiological factors such as the impedance of the electrode (Zaidel et al., 2010). To create homogenous PSD results for all recorded sites, the "relative" (normalized) power spectral density was calculated by dividing it by the total power of the signal between 0 and 3000 Hz. This relative, or normalized, power spectral density therefore estimates the spectral peak in relation to the other peaks in the spectrogram.
To compute coherence, the magnitude squared (MS-) coherence method (Kay, 1988;Miller and Sigvardt, 1998) was used. Welch's method was utilized, with a 1.5 s Hamming window (50% overlap), after removing the local window mean and with a spectral resolution of 1/3 Hz (same conventions as for PSD). Coherence values are limited (by definition) between 0 and 1. All coherence averages were therefore calculated in Fisher's transform domain (Miranda de Sa et al., 2009) and then reversed.
By definition, the removal of each window mean in the spectrum and the coherence analysis eliminate any power at 0 Hz (DC). We therefore start all the spectrum and coherence plot of this manuscript at 1 Hz.
A constant baseline level emerged in our coherence results (e.g., Figures 4A,B). This baseline probably resulted from the finite sampling of two "random noise" sources. To verify this, pairs of Gaussian random noise sources were simulated. The simulated data were subjected to the same filters and absolute operator as the real neuronal data and the same analysis tools. The magnitude of the coherence baseline dropped exponentially as time duration increased. Therefore, the baseline level seen in the STN coherence is most likely due to the finite (and relatively) short duration (mean = 23.7 s) of the recordings in human patients. The coherence functions were normalized by the subtraction of the average coherence of the randomly shuffled (10,000 times) pairs from the same STN domain. Note that the normalized coherence functions can therefore display negative coherence values.

SYNCHRONIZATION AND OSCILLATION STRENGTH (SYNC AND OSCIL SCORES)
Rosenberg and Halliday (Halliday et al., 1995(Halliday et al., , 2000Farmer et al., 1997;Farmer, 1998) proposed a very useful method to estimate coherence significance. However, this method employs a threshold confidence level, and does not offer a quantitative measure of synchronization strength. Therefore, a Z-score like method (effective Z-score: Z * ) was devised to determine the synchronization and oscillation strength. The Z-score of a given parameter is defined as the number of standard deviations above (or below) the mean. In this case the parameter was the maximum value (peak value) of the smoothed PSD or coherence (see below). However, instead of using the standard deviation of the entire frequency range, a tail standard deviation (σ tail ) was defined in the frequency range of 35 to 70 Hz. In this range, no coherence or power spectrum phenomena were observed in our dataset (Figures 1D, 2A-C, 3A,B,  4B). To smooth the coherence, a simple moving average (SMA) was calculated, with a window size of 23 samples (7.67 Hz), and a delta of one sample (i.e., the frequency resolution of 1/3 Hz). The synchronization strength or score was defined as is the maximum value of the coherence after smoothing with the moving average, and μ is the coherence mean. To find the frequency in which the spectrum or the coherence achieved maximal value, the arg-max(SMA(C(f ))) was calculated. The coherence (C(f )) maximal peaks were defined in the smoothed coherence function with minimal distances of 5 Hz between them. The search for the coherence peak was started at the lower frequency, and progressed to the largest value of the smoothed coherence function. All calculations (max, mean and arg-max) were performed in the frequency range between 1 and 70 Hz. Negative scores were found in a few cases, due to residual high power at the low (1-2 Hz) frequency range, and these were ignored.
To determine the oscillation strength, the same effective Zscore as for the synchronization was used and defined as the oscil score. The maximum value of the smoothed PSD (by a simple moving average, with window size of 23 samples, and delta of one sample, i.e., 1/3 Hz), and the tail standard deviation (σ tail ) in the frequency range of 35-70 Hz were calculated. The oscil score was defined as Z * = (MAX(SMA(PSD (f))) − μ/σ tail ). MAX(SMA(PSD (f))) is the maximum value of the PSD after smoothing with the moving average; μ is the PSD mean.
To explore the relationship between oscillation and synchronization a statistical measure of the oscillation strength of the two oscillatory sites was used. The average PSD * = (PSD1 + PSD2)/2 was calculated, where PSD1 and PSD2 were the power spectrum densities of each site in the neuronal pair, and the oscil score of PSD * was calculated. Additionally, other estimates of oscil scores of the two sites were calculated as: min(oscil 1 , oscil 2 ); max(oscil 1 , oscil 2 ); and as the geometric mean of the two scores where oscil 1 , oscil 2 were the oscil scores of each PSD site. The geometric mean was calculated as: oscil = sign(oscil1 * oscil2) * GeoMean(|oscil1|, |oscil2|), where sign was the sign operator of oscil 1 and oscil 2 product, in the case of negative values. Synchronization or oscillations were defined to be significant when the scores reached the Z * ≥ 2 (i.e., the coherence or the PSD peak value was higher than 2 SD of the mean values of these functions).

COHERENCE CONFIDENCE LEVEL
To assess the validation of sync score the confidence level (CL) of the coherence analysis (Halliday et al., 1995) was used. We divided the microelectrode records of duration R into L non-overlapped disjoint segments of duration S (R = L * S). The total spectrum was calculated using the average of the magnitude-squared (MS) of the discrete Fourier transform (periodogram), after removing the local mean in each segment S. Each segment contained S = 2 ∧ 16 samples with a frequency resolution of 0.7336 Hz. Only complete segments were analyzed; data points at the end of the record that did not make a complete segment were not included in the analysis. The procedures were implemented using Neurospec free MATLAB toolbox: http://www.neurospec.org. To obtain the approximate confidence interval for 95% and 99% from the data points, the level thresholds: CL 95 = 1 − 0.05 1/(L − 1) and CL 99 = 1 − 0.01 1/(L − 1) , respectively, were used. Figure 4B depicts examples of the relations between MS-coherence estimates (Z-scores) we used in the manuscript with coherence confidence levels of 95% and 99% respectively.

ASSESSING THE CAUSAL RELATIONS BETWEEN OSCILLATIONS AND SYNCHRONY
Spurious synchronization can arise from non-coupled oscillatory sites that oscillate in the same frequency bands (i.e., two atomic clocks might be synchronized due to their exact frequency although there is no physical coupling between them (Strogatz, 2003). To rule out this spurious oscillation-synchronization, the mean coherence of randomly shuffled pairs (10,000 times) was calculated for each category (all pairs, DLOR-DLOR, VMNR-VMNR, and DLOR-VMNR) of the STN. The shuffling was performed using the Mersenne Twister algorithm (Matsumoto and Nishimura, 1998) with a different seed number in each iteration. Figure 1 show an example of synchronous oscillatory activity as recorded by two electrodes inserted into the STN of a PD patient during DBS procedures. The raw analog data is shown in two time scales in Figures 1A,B. The power spectrums and the coherence function of this recording are shown in Figures 1C,D, respectively. One can easily observe the synchronous oscillations in the beta range (∼20 Hz) in this example.

SYNCHRONIZATION OCCURS ONLY BETWEEN DLOR PAIRS
To explore the properties of STN neuronal synchronization, STN spiking activity simultaneously recorded from two electrodes was analyzed (Figure 2). In total, 2390 multi-unit pairs along 72 STN trajectories (with >4 mm STN span in both electrodes) from 57 PD patients undergoing DBS surgery were included in the analysis. Previous physiological studies of the basal ganglia in the rodent (Mallet et al., 2008a,b) and primate (Bergman et al., 1994;Nini et al., 1995;Raz et al., 1996Raz et al., , 2000Raz et al., , 2001Goldberg et al., 2002) models of PD have indicated abnormal synchronicity of basal ganglia neurons as one of the major changes occurring in the network following dopamine depletion. Nevertheless, when the neuronal synchronization of simultaneously recorded STN sites (over the entire STN 2390 pairs) was measured, no distinguishable synchronization was found.
The STN can be spatially differentiated into sub-regions according to neural activity (Zaidel et al., 2010). Two areas could be robustly discriminated in our recording: the dorsolateral oscillatory region (DLOR, n = 1778 sites, Figure 3A) and the ventromedial non-oscillatory region (VMNR, n = 3002 sites, Figure 3B). Figures 3C,D show the distribution of the oscil scores in the DLOR and VMNR, respectively. As expected, significantly

FIGURE 2 | Synchronization found only between sites in the DLOR. (A)
Substantial coherence at beta frequency range exists between sites within the DLOR (from the STN entrance to the DLOR\VMNR border). However, the relative coherence is around zero in the beta range for pairs recording from the VMNR (B), or between sites of the DLOR and VMNR (C). HMM was used to automatically define the DLOR-VMNR border. The first row is calculated using the HMM defined border, and lower rows calculated with a progressively increasing gap from the HMM border. N is the number of pairs, and the shaded areas represent the standard error of the mean (SEM) of the relative coherence function for each population. (D) Distribution of oscillations scores along the normalized DLOR depth. For each trajectory the DLOR length was normalized for 0-entry, 1-end of DLOR.

Frontiers in Systems Neuroscience
www.frontiersin.org November 2013 | Volume 7 | Article 79 | 7 higher oscil score values were observed in the DLOR than in the VMNR. The division of the STN into the DLOR and VMNR domains enabled testing of the synchronization of STN pairs from the same and different regions. Significant synchronization, mainly in the frequency band of 8-30 Hz was found, but only between pairs in the DLOR itself (DLOR-DLOR, n = 615 pairs, Figure 2A, upper subplot). This synchronization was not observed in pairs of electrodes at DLOR and VMNR (n = 548 pairs, Figure 2B, upper subplot) or in the VMNR (n = 1227 pairs, Figure 2C, upper subplot). This finding is consistent with previous multiple electrode studies of the human Parkinsonian STN (Levy et al., 2000(Levy et al., , 2002aAmirnovin et al., 2004;Weinberger et al., 2006Weinberger et al., , 2009Alavi et al., 2013;Lourens et al., 2013) which reported coherence between STN oscillations in a small fraction of STN pairs. However, our findings indicated that the topographical location of the STN electrodes affected the probability of finding a correlation between STN sites, and coherence was only and robustly found between DLOR-DLOR multi-unit pairs.
Recent imaging studies (Lambert et al., 2012;Haynes and Haber, 2013) have clarified that the boundaries between the functional subdomains of the STN are fuzzy, and an overlap of motor and non-motor projections can be found in the transition areas between the STN domains. Therefore, the average coherence at the dorsolateral and ventromedial STN was tested with increasing gaps (0.5-2 mm) from the HMM borders. These results are shown in the lower five rows of Figure 2, and reveal a sharpening and increase of the average coherence peak in the STN DLOR when the gap is increased. Similarly, when the oscillation scores are calculated along the normalized depth of the DLOR, a gradual decrease in the oscillatory scores is observed as the DLOR lower border is approached (Figure 2D).
The above results were obtained by averaging over pairs recorded for different durations. The average coherence results (Figure 2) was further compared to the average coherence results of the same pairs with homogenous intervals (only the first 10 s of each recording was included, and recordings with durations shorter than 10s were excluded). Similar results (data not shown) were obtained.

SYNCHRONIZATION vs. OSCILLATIONS IN THE DLOR AREA
Next, correlation between the oscillations and synchronization in the STN was analyzed. The oscillation and synchronization strengths were calculated using the oscil and sync scores for each pair in the DLOR area. Figures 4A,B depict three examples of power spectrum and coherence function of STN activity with their relative oscil and sync scores, respectively. See also Figure 1 for an example of a simultaneous recording of two  Figure 5B depicts the scatter plot of the oscil and sync scores for all DLOR pairs. Different indicators for the oscillation strength of pairs of STN sites were used: the minimal and maximal oscil score, the arithmetic average of the PSDs, and the geometric mean of oscil scores. In all cases, the scatter plot of the sync score vs. the oscil score of the pairs within the DLOR area (n = 615) indicated a significant correlation (r > 0.24, p < 0.001) between the synchronization and the oscillations. Here (Figure 5B) we show only the data for the arithmetic mean of the oscil. scores.
The correlation between the oscillation and synchronization strength could imply that the synchronization pattern was dependent on the homogeneity of the neuronal oscillations within the DLOR. If the neural oscillations in different sites of the STN of a single patient have a very stable and equal frequency, the existence of synchronization may not be the result of physical coupling between the STN neurons. Therefore, the DLOR pairs of each trajectory were randomly shuffled and the synchronization between the shuffled (non-simultaneously recorded) pairs was re-quantified in each trajectory. After shuffling, the oscillations remained in the same frequency band, but the synchronization was no longer apparent. Figures 5A,C show the average coherence functions before and after shuffling of the DLOR-DLOR pairs FIGURE 5 | Synchronization is no longer apparent between non-simultaneously recorded (shuffled) STN DLOR pairs. (A) A significant average coherence between non-shuffled pairs within the STN DLOR. N = 615 represents the number of DLOR-DLOR pairs. (B) Scatter plot of synchronization (sync) and oscillation (oscil) scores in the STN DLOR reveals that the two measures are correlated. Each square represents the synchronization (Y -axis) vs. average (arithmetic mean) of the two oscillation scores of one of the 615 pairs within the DLOR. r is the Spearman correlation coefficient and p is the probability that r = 0 (no correlation between the scores). (C) Synchronization is no longer seen between non-simultaneously recorded (shuffled) STN DLOR-DLOR pairs. Inset: Schematic illustration of the shuffling procedure. The shuffling procedure was repeated 10,000 times for each pair. (n = 615), respectively. Figures 6A,B depict the average sync and oscil scores before and after shuffling in the STN DLOR and the VMNR. As expected, shuffling had no significant effect on the oscil score in either area (oscillation is a property of a single element and therefore should not be affected by the shuffling procedure). However, the average sync score of the DLOR pairs, but not the VMNR pairs declined significantly after the shuffling procedure ( Figure 6A). Finally, Figure 5D depicts the scatter plot of the sync and oscil scores of the shuffled pairs within the DLOR. The Spearman correlation between the sync score and oscill score dropped dramatically from r1 = 0.34 to r2 = 0.05 (p < 0.001 for the null assumption that r1=r2). The mean PSD estimate for the average oscil score (as in Figure 5B) was used for this analysis. Similar results were obtained for the other indicators of oscillation strength of the STN pairs.

COHERENCE IS MAINLY IN THE BETA FREQUENCY RANGE
Next, the frequency value where each spectrum ( Figure 7A) and coherence ( Figure 7C) reached its maximal value was calculated. In both cases, a bi-modal distribution was observed, with a dominance of tremor frequency (3-7 Hz) and beta (12-30 Hz) oscillations, for the auto-spectrums and the coherence functions, respectively. Figures 7B,D show the scatterplot of the maximal oscil and sync scores, respectively, as a function of their frequency. While the oscil scores had similar values in the beta and the tremor range, the values of the maximal sync scores were much higher in the beta than in the tremor range. These results are in line with our previous primate studies (Raz et al., 2000) that revealed mainly 5 Hz peaks in the auto-correlations vs. higher frequencies (10 Hz) in the cross-correlations functions of pallidal units and pairs recorded in the globus pallidus after MPTP treatment. However, we cannot rule out the possibility that the 10 Hz activity in this study is not tremor related and a harmonic feature (or n:m locking) of the tremor or of the neuronal oscillations at the tremor frequency.

LACK OF POSITIVE CORRELATION BETWEEN THE STN OSCIL AND SYNC SCORES vs. PD SYMPTOMS
Previous studies have suggested that STN oscillations and synchronization are correlated with tremor in PD patients (Levy et al., 2002a). This would indicate that the STN synchronized oscillations are driven by the tremor (which may be generated by an independent neuronal loop). The above findings of robust synchronization in the beta rather than in the tremor frequency range (Figure 7D) are not in line with this hypothesis. Nevertheless, we looked for correlations between the oscil and sync scores of our patients and their pre-operative (OFF medication) UPDRS scores. We did not find significant positive correlation between the average oscil score and sync scores of STN activity and the UPDRS scores of the tremor in the contra-lateral upper limb(s), all tremor (including axial) scores, and all UPDRS III motor scores. There is a trend for STN synchronized beta oscillations to be more robust in patients with less tremor. While these results might point to a correlation between STN beta oscillations and akinetic/rigid Parkinsonian symptoms (an issue that requires clarification in future studies with bigger sample of patients and with intra-operative clinical assessment), they definitely indicate that the STN beta synchronized oscillations are not a by-product of the PD tremor.

DISCUSSION
In this manuscript, synchronization within the human Parkinsonian subthalamic nucleus was investigated. No significant synchronization was found over the STN as a whole. After dividing the STN into two electro-physiologically distinct regions, the dorsolateral oscillatory region (DLOR) and the ventromedial non-oscillatory region (VMNR), significant synchronization in the beta range was observed, however, only within the DLOR. The strength of the DLOR synchronization was correlated with the strength of the oscillations of the multi-unit pairs. Nevertheless, shuffling between DLOR pairs abolished synchronization, suggesting that STN synchronization is an independent phenomenon and not a mere reflection of neuronal oscillations at similar frequencies.
Previous studies have shown significant spatial overlap between the DLOR and the STN sensorimotor area (Rodriguez-Oroz et al., 2001;Zaidel et al., 2010). The finding that the STN VMNR (considered to be part of the limbic and associative basal ganglia network) remains unsynchronized is consistent with the predominantly motor nature of PD. However, the (normal) lack of synchronization in the STN VMNR may be due to a selection bias of our DBS patients. Since conventional inclusion criteria were used to select candidates for DBS, patients were usually severely motor-impaired and had few of the non-motor features of the disease. Furthermore, the DLOR may reflect the pathological area of the STN which progressively invades the limbic domains of the STN as the disease advances. Finally, our results are in line with a fuzzy rather than a sharp boundary between the STN sub-domains (Lambert et al., 2012;Haynes and Haber, 2013).

THE STN SPIKING POPULATION ACTIVITY IS SYNCHRONIZED
In this study population spiking (multi-unit) activity was used as a measure of the spiking activity of the STN rather than the more classical parameter of single unit activity (Perkel et al., 1967;Abeles, 1982;Lemon, 1984;Eggermont, 1990). This was primarily for practical reasons. The goal of physiological recording in the operating room (OR) is to enable better identification of the borders of the subthalamic nucleus and its sub-regions. The electrode is therefore advanced in 100 μm steps rather than 2-5 μm steps as is customary in the research laboratory setup. The sampling duration at each step is also limited (Shamir et al., 2012) and the OR conditions often do not allow stable recordings (as compared to 30-90 min stable recordings in a research setting). On the other hand, the cross-correlation of composite spike trains derived from several un-discriminated cells recorded on a single electrode (multi-unit activity) enhances the sensitivity of correlation methods. First, the higher discharge rate of multi vs. single unit recording reduces the asymmetric sensitivity of correlation methods to excitation vs. inhibition (Aertsen and Gerstein, 1985). Second, multi-unit cross-correlation can be a more sensitive detector of a neuronal relationship than single-unit cross-correlation (Bedenbaugh and Gerstein, 1997 reasons. Furthermore, the use of signals recorded by two different electrodes in this study reveal the long range (2 mm) synchronization of STN DLOR. It is hoped that future studies of STN units using objective metrics for quantification of the quality of the unit isolation (Joshua et al., 2007;Hill et al., 2011) will shed more light on synchronization in the STN and other basal ganglia structures of human patients.

SYNCHRONIZATION ONLY OCCURS BETWEEN DLOR PAIRS IN THE STN
Early studies described neuronal synchronization in the STN as an epiphenomenon found mainly in patients presenting with tremor (Levy et al., 2000(Levy et al., , 2002a. More recent studies (Hanson et al., 2012;Alavi et al., 2013;Lourens et al., 2013) have reported that synchronization can be found between some but not all STN pairs. On the other hand, beta-band LFP and spike oscillations have been described as a consistent feature of human PD in the dopamine depleted state (Brown and Williams, 2005;Foffani et al., 2005;Little et al., 2012). Moreover, many studies have documented the consistency of beta-band oscillations in both the spatial and temporal domains (Bronte-Stewart et al., 2009;de-Solages et al., 2010;Zaidel et al., 2010;Abosch et al., 2012;Little et al., 2012). In this study, synchronization within the Parkinsonian STN DLOR was indeed found to correlate with oscillations. However, the shuffling procedure revealed that STN synchronization was not due to independent oscillators with a similar oscillation frequency (Strogatz, 2003). If this had been the case, a significant synchronization should also have been observed between the shuffled (non-simultaneously recorded pairs of the same patient) DLOR-DLOR pairs. Thus, the synchronization of the simultaneously recorded STN pairs probably reflects the increased coupling between these neurons in the dopamine depleted state of Parkinson's disease. This increased coupling is probably due to the increased efficacy of the common inputs to the STN cells, either from the cortex (Nambu, 2004;Kita and Kita, 2012) or from the external segment of the globus pallidus (Plenz and Kitai, 1999;Tachibana et al., 2011). However, at this stage the possibility of increased coupling by lateral connectivity within the STN cannot be ruled out (Parent et al., 2000;Parent and Parent, 2007). The finding that most of the energy of the STN synchronous oscillations is in the beta range suggest that these oscillations are not generated by feedback of the peripheral tremor. It is interesting to note that synchronous oscillations in the basal ganglia of MPTP treated primates are mainly found in the 10 Hz domain, where human oscillations span the full beta range (12)(13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30). Future studies should reveal if this is due to species difference, or due to differences between the MPTP model and human idiopathic Parkinson's disease.

CONCLUDING NOTES
In this study we show that the STN domain most affected by PD dopamine depletion (the DLOR, probably the STN motor domain) exhibited both oscillations and synchronization. This suggests that synchronization reflects an additional property of the Parkinsonian STN. Previous studies in the basal ganglia of MPTP treated primates have demonstrated that synchronization can be completely independent of oscillatory activity (Heimer et al., 2002). The previous and the current findings can serve to recast the relationship between oscillations and synchronization in the Parkinsonian basal ganglia (Raz et al., 2000;Amirnovin et al., 2004;Moran et al., 2008). In addition to changes in discharge rate and pattern, synchronization within the STN may be another pathophysiological marker of Parkinson's disease. The potential consequences of synchronization (as opposed to other attributes like rate and pattern change) are probably mainly due to reduced information capacity of the basal ganglia neurons. However, the different pathological changes in the parkinsonian basal ganglia are probably not mutually exclusive. Synchronized oscillations have stronger effects than less synchronized oscillations and completely unsynchronized oscillations might have no effect on target neurons. Furthermore, future studies toward adaptive DBS (Rosin et al., 2011) should investigate which of the pathophysiological changes in the STN activity might be used as the optimal trigger for closed loop DBS.

ACKNOWLEDGMENTS
This work was supported partially by the Vorst Family Foundation for Parkinson Research, by the Simone and Bernard Guttman chair of Brain Research, and the generous support of the Rosetrees and Dekker foundations (to Hagai Bergman) as well as by the ELSC post-doctorate fellowships to Shay Moshel and Reuben R. Shamir and the PATH and Bloom foundations of London to Zvi Israel. We would like to thank Prof. Y. Ritov for statistical advice and Dr. S. Freeman and E. Singer for help in scientific editing.

AUTHOR CONTRIBUTIONS
Shay Moshel and Reuben R. Shamir claim for equal contribution. Shay Moshel did the data analysis, and wrote the manuscript together with Hagai Bergman and Zvi Israel. Reuben R. Shamir handled the database and collected part of the data. AR initiated the study and helped in the collection of data. Hagai Bergman, Renana Eitan, Fernando R. de Noriega, and Zvi Israel collected the data. All authors discussed the results, reviewed the manuscript, and made their comments.