Unsupervised clustering reveals spatially varying single neuronal firing patterns in the subthalamic nucleus of patients with Parkinson's disease.

INTRODUCTION: Subthalamic nucleus (STN) is an effective target for deep brain stimulation (DBS) to reduce the motor symptoms of Parkinson's disease (PD). It is important to identify firing patterns within the structure for a better understanding of the electro-pathophysiology of the disease. Using recently established metrics, our study aims to autonomously identify the discharge patterns of individual cells and examine their spatial distribution within the STN.
METHODS: We recorded single unit activity (SUA) from 12 awake PD patients undergoing a standard clinical DBS surgery. Three extracted features from raw SUA (local variation, bursting index and prominence of peak) were used with k-means clustering to achieve the aforementioned unsupervised grouping of firing patterns.
RESULTS: 279 neurons were isolated and four distinct firing patterns were identified across patients: tonic (11%), irregular (55%), periodic (9%) and non-periodic bursts (25%). The mean firing rates for irregular discharges were significantly lower (p < 0.05) than the rest. Tonic firings were significantly ventral (p < 0.05) while periodic (p < 0.05) and non-periodic (p < 0.01) bursts were dorsal. The percentage of periodically bursting neurons in dorsal region and entire STN were significantly correlated with off state UPDRS tremor scores (r = 0.51, p = 0.04) and improvement in bradykinesia and rigidity (r = 0.57, p = 0.02) respectively.
CONCLUSION: Strengthening the application of unsupervised clustering for firing patterns of individual cells, this study shows a unique spatial affinity of tonic activity towards the ventral and bursting activity towards the dorsal region of STN in PD patients. This spatial preference, together with the correlation of clinical scores, can provide a clue towards understanding Parkinsonian symptom generation.

Introduction

The subthalamic nucleus (STN) is an established target in deep brain stimulation (DBS) to control the motor symptoms of Parkinson's disease (PD). It has been suggested that firing pattern changes in single unit activity (SUA) within the basal ganglia could be a key in understanding the pathophysiology of PD [1] and can assist in electrophysiological mapping of the structure [2]. Studies involving 1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine (MPTP) treated non-human primates report an increase in oscillatory neuronal activity and instantaneous firing rates [3] along with an increase in bursting activity [4]. Based on these observations, Wichmann and DeLong suggested that both the change in firing rate and firing patterns are relevant to the pathophysiology of PD and investigating them have the potential to unravel the underlying neural network inducing PD manifestations [5]. In this regard, various groups have reported different types of neuronal firing patterns in the STN of Parkinsonian patients [2,3,6,7]. In particular, oscillatory firings were separated from non-oscillating ones in [7], whereas others have focused on identification of tonic (regular), irregular and oscillatory firing patterns [6]. The majority of these studies characterized firing profiles in a binary fashion. The metrics used to separate these profiles are modestly unclear and, in several instances, chosen arbitrarily based on visual inspection such as inter-spike interval (ISI) histograms or binary spike trains. Notwithstanding, in total, three main firing patterns have been commonly described: tonic, irregular and oscillatory [2,3,6,7]. Irregular firings are most common in proportion, followed by tonic and then oscillatory [2,6]. Furthermore, studies also investigated the spatial distribution of oscillatory firing patters [7] as the dorsolateral STN is known to harbor the motor circuitry of STN, whereas the ventromedial STN is associated with limbic circuitry [6].

Grouping of firing patterns into different categories manually is difficult as it is both time consuming and, in some cases, hard to achieve using the limited information provided by an auditory/visual inspection. By using recently established metrics [8,9], our study aims to automatically characterize the proportion and spatial distribution of four firing patterns: tonic (T), irregular (I), periodic bursts (PB) and non-periodic bursts (NPB) within the STN of PD patients revealing distinct firing profiles in the territories of STN.

Materials & methods

Patients and recordings

The study was approved by the Institutional Review Boards at the University of Houston and Baylor College of Medicine. Intraoperative SUA was recorded from 12 awake PD patients (8 Male, 4 Female, Age: 54 ± 11.14 years) under local anesthesia as part of the standard clinical procedure. Eight of the 12 patients underwent bilateral DBS implantation, and 4 had unilateral implants constituting 20 STN recordings in total. Patients were asked to stop medication at least 12 h prior to the surgery. Three microelectrodes (NeuroProbe; AlphaOmega, Israel) were initially placed at least 15 mm above the stereotactic target and advanced towards the target depth with 0.5 mm steps. SUA was recorded from three micro-electrode trajectories using a bio-amplifier (Grapevine Neural Interface and Processor; Ripple, UT) at 30 KHz. During the awake surgery, real time SUA traces were displayed on a data-scope to help with the localization of the target structure. The entry and exit of the STN was determined by a clinical neurophysiologist by listening to and visually observing the firing patterns of neurons. The dorsal borders of the STN in each patient were defined by inspection of SUA showing an increase in background activity and discharge rates. The dorsal border was taken as reference and marked as 0 mm with negative depth values ventral to it. Spatially, 0 mm to −2.5 mm was considered as the dorsal half of STN as previous studies show dominance of oscillatory activity in this area [7]. Other firings captured above or below the electrophysiologically defined borders of the STN were not included in our report. The de-identified data were transferred to a computer for offline processing.

Data analysis

The neural data was analyzed offline in Matlab 2017b (Mathworks, MA) using custom scripts and a publicly available spike sorting toolbox [10]. Since the neurons are observed to display a non-stationary firing behavior over time [2], we analyzed and clustered the raw segments from all 3 tracks in 5 s intervals with a 2.5 s overlap. This procedure helped us in capturing the neurons' non-stationary behavior and label each segment separately based on its firing characteristic. Consequently, rather than assigning each neuron to a single category such as tonic or irregular, we clustered the segments and assigned them into categories.

Following the high pass filtering of raw SUA data over 300 Hz, spike detection and sorting were performed automatically by using WaveClus [10], an open-source Matlab toolbox. Segments with a firing frequency of <10 Hz and duration <3 s were discarded. Three features: (i) local variation [8]; (LV), (ii) bursting index [9]; (BI) and (iii) prominence of peak (P) were used together with k-means [11] to group the spike trains. To analyze the spatial distribution of discovered clusters statistically, Wilcoxon sign rank test was performed within each firing pattern to check if the pattern was significantly dorsal or ventral and to assess the difference in firing rates between these patterns as well as the intensity of bursting (BI) and regularity (LV) in dorsal and ventral regions. Correlation with Unified Parkinson's Disease Rating Scale (UPDRS) was done using Spearman rank correlation due to non-normality of the data. Three patients were excluded from the correlation analysis since they were external referrals and access to their detailed medical records and itemized UPDRS scores were limited (see supplementary Table S1).

Local variation (LV)

LV was used as a measure of regularity in firings and defined as:where Ii and Ii+1 are the ith and i + 1st ISI. LV helped eliminate the issues caused by variable firing rates and non-stationarity [8]. LV approaches 0 for a tonic firing neuron and 1 for a highly irregularly firing neuron.

Bursting index (BI)

Bursts were detected using a modified version of Max Interval Method [9]. For burst initiation, ISI below 10 ms and a minimum of three consecutive spikes were required along with an intra-burst pause period of at least 20 ms or a sum of 30 ms between the next two ISIs. The bursting index is accepted as the ratio of the number spikes bursting to the number of spikes with ISI >10 ms.

Prominence of peak (P)

Upon initial visual investigation of the raw SUA data we observed that bursting can be both periodic and non-periodic. In order to distinguish between periodic and non-periodic bursts, a new feature called prominence of peak was computed from the spike train power spectrum. The digital spike train was convolved with a 10 ms Gaussian window and its power spectral density (PSD) then calculated using Welch periodogram [12]. The peak value of the PSD was detected between 3 and 40 Hz and its prominence was calculated using:where P is the prominence, P is the peak value of PSD and is the mean of the minimum values around this peak within an 8 Hz bandwidth. A regularization parameter (r = 10−3) was used to eliminate small peaks that might appear prominent due to a small background. Eq. (3) was used to limit the prominence between 0 and 1 where values approaching 0 indicate strong oscillatory activity.

Clustering

Clustering of firings into four categories was done using the k-means method [11] in a 3D feature space. The features were normalized before being subjected to the clustering algorithm. We varied the number of clusters from 1 to 8 and computed the sum of within-cluster distances for each index. When a relatively prominent elbow was observed in the graph, we used that index to decide on the total number of clusters in the data.

Results

We analyzed 6705 s of data and successfully isolated 279 neurons (23 ± 7 neurons per patient). After spike detection and sorting, neurons were grouped into firing patterns and mapped onto the STN to investigate their spatial distribution. Fig. 1A shows the border and depth definitions as well as exemplary raw SUA recordings from various depths. Fig. 1B displays a representative trace for each irregular, tonic, periodic and non-periodic bursting firing pattern we studied. The clustering performed using the three features provided a prominent elbow around four for the total sum of intra-cluster spacing suggesting four distinct firing patterns (Fig. 1C). To provide an example of their content, we show 3 raw traces from each cluster as well. We observed that the largest cluster corresponding to 55% of the neuronal spiking displayed irregular patterns, where the remaining 11% displayed a tonic pattern while other 34% were bursting with 25% non-periodically and 9% periodically (Fig. 1D). The mean firing rates for irregular firings (33.48 ± 17.7 Hz) were significantly lower (p < 0.05) than the tonic firings (74.19 ± 33 Hz), periodic (82.35 ± 43.46 Hz) and non-periodic bursts (57.57 ± 28.66 Hz). Visual observation of the spatial distribution of these firing profiles (Fig. 1E) revealed that the irregular firings were uniformly distributed throughout the STN. Interestingly, burst type firings populated the dorsal STN whereas tonic ones tended towards the ventral end. We observed that the amount of tonic activity rapidly increased 3 mm below the dorsal border of STN and extended towards the ventral region. Statistical analyses concluded that tonic firings were significantly ventral (p < 0.05) while periodic (p < 0.05) and non-periodic bursts (p < 0.01) were significantly dorsal. However, it should be noted that bursting patterns can still be observed within the ventral territory of the STN. Similarly, our findings show that both bursting (BI: 0.86 ± 0.17 vs 0.79 ± 0.2) and regularity (LV: 0.31 ± 0.06 vs 0.21 ± 0.08) were significantly different in the two regions of STN (p < 0.01) with bursting being more intense in the dorsal region and tonic discharges being more regular in the ventral region. Finally, as done in [7] we correlated the off drug UPDRS part III scores and their improvement with medication with firing characteristics and found that the percentage of dorsally located periodic bursting neurons was significantly correlated with tremor (r = 0.51, p = 0.044). The percentage of periodic bursting neurons in the entire STN was significantly correlated with the improvement in bradykinesia and rigidity (r = 0.57, p = 0.02). Other UPDRS sub-scores did not show any significant correlation (p > 0.05).

Fig. 1

(A) Representation of STN and raw SUA recordings from various depths (0 mm is the reference dorsal border). (B) Exemplary raw SUA data for four firing patterns (I = Irregular, T = Tonic, PB = Periodic burst, NPB = Non-periodic burst). (C) Clustering result over 12 subjects (1371 segments) (LV = Local Variation, BI = Bursting Index). Inset - Inter cluster distance plots for 1 to 8 clusters (x-axis: number of clusters, y-axis: sum of inter-cluster distance) and three examples of raw SUA traces from each cluster. (D) Proportion of firing patterns distributed across 12 subjects. (E) Spatial distribution of firing patterns across STN. The shaded areas indicate dorsal (0 to −2.5 mm) and ventral (−3.5 to 6 mm) regions of STN. (* Indicates significant difference between shaded regions) (* = p < 0.05, ** = p < 0.01).

Discussion

In this study, using an unsupervised approach utilizing recently established metrics, we identified and investigated the spatial distribution of neurons within the STN of PD patients based on their firing patterns. In addition to commonly reported irregular, tonic and bursting patterns [2,6], we further divided the bursting category into oscillatory and non-oscillatory bursts. While most studies compare the spatial distribution of oscillatory units against non-oscillatory units [7], we expanded the investigation by looking at the spatial distribution of each of the four individual patterns separately. An interesting and novel finding was the occurrence of tonic firings towards the ventral region of STN. While the presence of oscillatory neurons in the dorsal region has already been established [7], our study shows the occurrence of burst type of findings (be it periodic or non-periodic) in the dorsal region as well. Our observations of irregular firings dominating the set (55%), followed by non-periodic bursts (25%), tonic (11%) and periodic burst firings (9%) were comparable to previous findings [2,3]. The high fraction of bursting in our findings is also in accordance with the previously studied MPTP treated animal models [3,4] where the main difference between healthy and Parkinsonian state was the emergence of bursting activity and increased instantaneous firing rates. Similar to the study conducted by Weinberger et al. [7], we found a significant correlation between the percentage of periodically bursting dorsal neurons with the off state tremor score and overall percentage of periodic neurons with improvement in bradykinesia. Given the somatotopic organization of STN [6], our observation of the predominance of bursting activity in dorsal STN and tonic activity in ventral region could be indicative of variable disease manifestations such as motor and cognitive/limbic features of PD. The study of pattern changes as well as their distribution within the structure thus could provide a better understanding of the more complex mechanisms underlying the disease.

We presented an automated approach that can isolate clusters of neurons based on their firing profiles in the territories of STN. In practice 83% of centers around the world are using MER for electrophysiological mapping to refine the final position of the chronic electrode [13]. Intraoperative decision making involves the conversion of single-unit neuronal activity recorded at the microelectrode tip, into audio and visual signals. This process of mapping is time consuming experience-based and depends critically on the ability of neurosurgeons or electrophysiologists to recognize a variety of cues in the recorded neural activity. Multiple microelectrodes isolating neurons with distinct firing patterns makes the interpretation of auditory cues even more challenging. We propose that the presented automated approach can be used to investigate the firing profiles of a large number of neurons recorded simultaneously during DBS surgery. We postulate that it may have utility in identifying the most favorable track for macroelectrode implantation, and contact localization (i.e. depth), while also providing additional insight regarding the relationship of disease symptoms to different territories of the STN. Since the bursting firings were biased towards the dorsal region, the track with larger number of bursting neurons can be utilized for chronic DBS electrode placement. Nevertheless, further research will be required to better define how this technique can be applied to improve targeting of the dorsolateral STN where the majority of patients will experience clinical improvement.

It must be noted that using existing clinical electrodes one can record only 1–2 neuron at best per track at each depth thus limiting the number of neurons isolated within the territories of STN of each patient. We anticipate that, an increased patient population as well as future clinical electrodes designs with multiple SUA recording contacts can improve the sampling capacity and will likely provide advantages in establishing a bridge between single neuronal activity and clinical observations.

Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Declaration of competing interest

None.