Pre-print manuscript from:
Spatio-temporal dynamics of cortical drive to human subthalamic nucleus neurons in Parkinson's disease
(Sharott et al., 2018)

Summary .
Pathological synchronisation of beta frequency (12-35Hz) oscillations between the subthalamic nucleus (STN) and cerebral cortex is thought to contribute to motor impairment in Parkinson’s disease (PD). For this cortico-subthalamic oscillatory drive to be mechanistically important, it must influence the firing of STN neurons and, consequently, their downstream targets. Here, we examined the dynamics of synchronisation between STN LFPs and units with multiple cortical areas, measured using frontal ECoG, midline EEG and lateral EEG, during rest and movement. STN neurons lagged cortical signals recorded over midline (over premotor cortices) and frontal (over prefrontal cortices) with stable time delays, consistent with strong corticosubthalamic drive, and many neurons maintained these dynamics during movement. In contrast, most STN neurons desynchronised from lateral EEG signals (over primary motor cortices) during movement and those that did not had altered phase relations to the cortical signals. The strength of synchronisation between STN units and midline EEG in the high beta range (25 - 35 Hz) correlated positively with the severity of akinetic-rigid motor symptoms across patients. Together, these results suggest that sustained synchronisation of STN neurons to premotor-cortical beta oscillations play an important role in disrupting the normal coding of movement in PD.
Introduction.
Synchronised oscillations in the cortico-basal ganglia network are a prominent feature of patients with Parkinson's disease (PD) and its animal models (Brown et al., 2001; Sharott et al., 2005; Tachibana et al., 2011). Oscillatory activity in the beta-range (broadly defined here as any frequency between 12-35 Hz) in the STN of PD patients is coupled to that in the cortex and basal ganglia output nuclei (Brown et al., 2001; Williams et al., 2002). The power and synchrony of these signals decreases following therapeutically effective dopamine replacement and deep brain stimulation (Marsden et al., 2001; Williams et al., 2002; Kühn et al., 2008). Consequently, cortico-subthalamic beta oscillations, and those in the wider basal ganglia network, have been proposed as a putative pathophysiological mechanism in PD (Williams et al., 2002; Eusebio and Brown, 2009).
Elucidating the mechanism by which pathological beta oscillations are generated and propagated through different parts of the network is crucial to advancing this hypothesis (Brown, 2007; Stein and Bar-Gad, 2012). Frontal/midline EEG signals (assumed to reflect the supplementary motor area) lead STN local field potentials (LFPs) with a time delay of 20- 30 ms (Williams et al., 2002; Fogelson et al., 2006). In contrast, cortical beta oscillations recorded using ECoG electrodes over primary motor cortex led STN phase-locked spiking, revealed using spike trigged averages, by around 100 ms (Shimamoto et al., 2013). Due to differences in recording and analysis methods in these studies, it is not clear to what extent the dynamics of synchrony differ between single units and LFPs, and/or between different cortical areas.
During movement, cortico-subthalamic LFP beta desynchronisation correlates with motor performance (Kühn et al., 2004; Tan et al., 2015; Tan et al., 2016). Such task related changes in oscillation and synchronisation could reflect a role for these activities in healthy movement (Tan et al., 2016). In contrast, traditional models of basal ganglia function, based on recordings of single neurons in the STN and other BG nuclei, focus on widespread changes of the firing rate of individual neurons in coding movement parameters (Alexander and Crutcher, 1990; DeLong, 1990; Mink, 1996). It is unclear if, and how, such movement related changes in STN unit firing rate and pattern relate to beta synchronisation. Specifically, to what extent do individual STN neurons uncouple from cortical beta oscillations during movement-related firing rate changes? If STN neurons do not completely uncouple from cortical activity, does the angle of phase locking change in comparison to resting activity? Such questions have assumed greater significance due to the development of closed-loop deep brain stimulation, where the time delays between cortical oscillations and basal ganglia stimulation are an important determinant in the ability of these approaches to reduce motor symptoms (Rosin et al., 2011).
Here we demonstrate that the spikes of individual STN neurons are locked to cortical beta oscillations with constant time delays that differ between frontal/midline and lateral cortical areas. During periods of self-paced movement, many units only partially desynchronise from frontal and midline cortical beta oscillations and maintain the same phase locking dynamics that occur during rest. In addition, the temporal and spatial extent to which STN units are phase locked to midline cortical beta oscillations is highly correlated with the patients’ akinetic/rigid symptoms. Together, these results show that the dynamics of beta synchronisation between STN neurons and premotor-cortical beta oscillations are highly stable and could contribute to motor impairment in PD.

Methods.
All statistical values reported in the methods section are given as mean ± standard deviation (SD) unless noted otherwise.
Patient details and clinical scores .
This study was conducted in agreement with the Code of Ethics of the World Medical Association (declaration of Helsinki, 1967) and was approved by the local ethics committee. All patients gave their written informed consent to participate in this study. We studied 12 patients (7 female, 5 male; age: 65 ± 6 years) suffering from advanced idiopathic PD (Hoehn & Yahr score: 3.5 ± 0.7, (Hoehn and Yahr, 1967) with a disease duration of 16 ± 5.8 years. Patients underwent bilateral microelectrode-guided implantation of DBS electrodes in the STN. Preoperatively, all patients had significant improvement of the Unified Parkinson’s Disease Rating Scale motor section score (part III; (Fahn et al., 1987) following intake of levodopa (symptom improvement: 53% ± 14). Furthermore, all patients were classified as cognitively intact (based on their performance on the Mattis Dementia Rating Scale (Mattis, 1988) and fulfilled other inclusion criteria for STN-DBS, such as no structural alterations on magnetic resonance imaging (MRI), and no concomitant severe medical comorbidities. Further clinical details are summarised in Supplemental Table 1. The patient ́s motor UPDRS scores were assessed one week (5 ± 3 days) before the operation, both OFF and ON medication (motor score in Dopa-OFF condition: 39 ± 13 vs. Dopa-ON: 19 ± 9; p < 0.001, paired t-test).


Cortical EEG and ECoG signals .
In parallel with the deep brain signals, a 32 channel system (AlphaLab, Alpha Omega Inc., Nazareth, Israel) was used to amplify and record electroencephalographic (EEG, amplification x4.000; bandpass, 0-400 Hz) and electromyographic (EMG; amplification x2.000; bandpass, 5-1000 Hz) signals at a sampling rate of 3005 Hz during the microrecording and stimulation periods (Moll et al., 2015). EEG was recorded from 4 scalp electrodes (Ag/AgCl cup electrodes filled with conductive gel; Nicolet Biomedical, Madison, WI, USA) placed approximately at positions Fz, Cz, C3 and C4 according to the international 10-20 system for the placement of EEG electrodes, with the left earlobe as a common reference. In 7 patients, we additionally recorded an electrocorticogram (ECoG) from the dura above the craniotomy, approximately over the dorsolateral prefrontal cortex and referenced to the Cz channel (Fig. 1A-B). All EEG channels were also re-referenced to the Cz electrode to give three bipolar EEG channels which will be referred to from here on as midline (for the Fz/Cz derivation), and contralateral or ipsilateral (for the Cz to C3/C4 derivation in respect to the side of the depth micro- or macroelectrode being analysed). In one case, the Fz position was not recorded, and has therefore only been included in analyses of ECoG.
Data selection and processing .
All data were collected as part of our routine neuronavigation process which aims to precisely identify the surgical target (STN). Patients were only included in this study if a minimum of 10 STN unit recordings, each lasting more than 45 s, were recorded from one or both hemispheres (Sharott et al., 2014). It is important to note that recordings were only included if they were made during sustained periods where patients were awake and highly co-operative. In total, 269 STN unit recordings were used in this study. In some recordings (96/269 units recorded in 9/12 patients), patients were engaged in a simple, brief movement task (flexion/extension of wrist or ankle) as part of the routine mapping procedure.
Spike detection was performed offline using a voltage threshold method, the threshold value of which was set sufficiently high relative to the noise level to avoid false-positives (typically, threshold values > 4 SD of the background level were used). When possible, single unit activities (SUAs) were then separated by manual cluster selection in 3D feature space on the basis of several waveform parameters including principal components, signal energy, peak time and the presence of a central trough in the autocorrelogram (Offline-Sorter, Plexon Inc., Dallas, TX, USA). Over half of the unit activities were characterized as SUAs (n = 184). Given the high density of neurons in the STN, we chose to be conservative in spike sorting and in the classification of single units. Many recordings classified as multiunit activities (MUAs) may, therefore, have comprised mostly of one or two STN units and were clearly distinct from background spiking activity utilised in some previous studies (Moran et al., 2008; Zaidel et al., 2010). In all cases, the initial 5-15 s of the unit recording was discarded as were any portions of injury discharge at the end of the recording.

Signal processing .
The firing rate of each neuron was calculated in the most stable part of the spike train. STN units are autonomously active (Do and Bean, 2003) and their tonic firing does not generally fall below 10 spikes/s in well isolated recordings (Magill et al., 2000; Mallet et al., 2008b). The firing rate of each neuron was calculated using the longest part of the spike train where the firing frequency was above 10 spikes/s in every 2 s bins for more than 30 s. By doing this, we attempted to ensure that mean firing rate was not affected by recording instability. The mean rate of single units was 33.7 spikes/s, in agreement with previous studies with well-isolated units in parkinsonian animals and humans (Bergman et al., 1994; Weinberger et al., 2006; Mallet et al., 2008b; Sharott et al., 2014).
Spectral parameters for both time-series and point-processes were evaluated using fast Fourier transform (FFT) as described in Halliday et al. (1995). A Hanning window filter was used for all spectral analyses and spectra were estimated by averaging across these discrete sections (Halliday et al., 1995). Spectral analyses and related correlations were computed with a combined population of single- and multi-units and were calculated using the MATLAB toolbox Neurospec (www.neurospec.org). Coherence and phase analysis were used to evaluate the strength and temporal delay of coupling between the EEG/ECoG and STN unit/LFP channels respectively. The cross- and autospectra from two time series can be combined to give the coherence: a measure of the degree to which one can linearly predict change in one signal given a change in another signal (Brillinger, 1981; Halliday et al., 1995; Rosenberg et al., 1989). Being a normalized measure, coherence is bounded from 0 to 1, with a value of 0 indicating non-linearly related signals and a value of 1 signifying two linearly identical signals. Significance was evaluated using confidence limits based on the number of segments used and was independent of frequency (Halliday et al., 1995). The variance of the coherence at each frequency was normalized by transforming its square root using the Fisher transform (Brown et al., 2001; Sharott et al., 2005).
When coherence between time series/point processes is significant, the phase spectra, defined as the argument of the cross-spectrum, can be used to calculate the time delay between two time series over a given frequency range (Halliday et al., 1995; Fogelson et al., 2006; Magill et al., 2006). Using a series of contiguous frequency bins, rather than a single frequency point to evaluate the phase delay is advantageous, as the latter is ambiguous (Halliday et al., 1995; Fogelson et al., 2005). Over frequency ranges with significant coherence (> 6 contiguous bins over the 95% confidence interval), the time lag/lead between two signals was calculated from the slope of the phase estimate, having fitted a line by linear regression to establish the slope as significantly different from zero (Halliday et al, 1995; Fogelson et al, 2005). A significant positive or negative gradient indicates a lag or lead by the input channel respectively relative to the reference channel, at a latency given by the slope of the gradient divided by 2 (Halliday et al., 1995). The time delay was calculated when the phase gradient (corresponding to the significantly coherence bins) had p-value of less than 0.05 and an R2 value of more than 0.5.

Phase locking analysis .
To further investigate how the activity of STN units varied in time with respect to ongoing cortical activity, we analysed the instantaneous phase relationships between STN spike times and cortical oscillations in narrow frequency bands. EEG, ECoG and LFP signals were first filtered, using a neutral-phase bandpass filter (Butterworth filter, 2nd or 3rd order) in 13, 5 Hz wide, frequency bands from 5 to 40 Hz with a 2.5 Hz overlap. Subsequently, the instantaneous phase and power of the ECoG in these frequency bands were separately calculated from the analytic signal obtained via the Hilbert transform (Lachaux et al., 1999). In this formalism, peaks in the oscillation correspond to a phase of 0° and troughs to a phase of 180°. Circular phase plots and circular statistical measures were calculated using the instantaneous phase values for each spike. Descriptive and inferential circular statistics were then calculated using the CircStat toolbox (Berens, 2009) for MATLAB. Units were first tested for significantly phase-locked firing (defined as having p < 0.05 in Rayleigh’s uniformity test). The null hypothesis for Rayleigh’s test was that the spike data were distributed in a uniform manner. For each of the neurons that fulfilled these criteria, the mean phase angle was calculated. The mean resultant vector length (referred to hereafter as simply ‘vector length’) of the phase distribution, bound between zero and one (the closer to one, the more concentrated the angles), was used to quantify the level of phase locking around the mean angle for individual neurons (computed using the angles of each spike) and for populations of neurons (computed using the mean angle for each neuron).

Analysis of rest and movement periods .
Movement periods, consisting of continuous flexion/extension of the wrist contralateral to the recorded STN, were isolated using EMG recordings from the flexor and/or extensor muscles. All movement periods were examined by eye, and those with 2 or more sustained rest (11.2 ± 6.5 s) and movement periods (17.8 ± 8.2 s), with a minimum of 4s, were defined using a threshold set by visual inspection. 42 STN units from 5 patients, of which 30 were SUAs, where recorded during such rest/movement periods. For each unit, the firing rate was calculated across movement and rest periods and converted to percentage change during movement. As we were interested in the relationship between firing rate and oscillation, for this analysis we normalised the mean resultant vector length against that of 100 ISI-shuffled surrogate spike trains to ensure the number of spikes could not confound the phase locking measure (Rivlin-Etzion et al., 2006; Sharott et al., 2009; Sharott et al., 2014). The vector length was therefore expressed as a z-score (the number of standard deviations of the vector length of the real spike train compared to this surrogate distribution). In order to compare with the % firing rate change, the phase-locking z-score for movement was subtracted from the rest z-score to give a single desynchronisation value. To evaluate the power of the cortical signals during these rest and movement periods, the amplitude envelope from the Hilbert-transformed beta-filtered data (at the frequency at which there was a phase-locked unit) for the entire rest and movement period was normalised to a z-score. The mean amplitude of the rest and movement periods was then calculated from this normalised data, allowing comparison of signal amplitude in rest and movement periods across recordings from different patients.
Correlation of phase locking and clinical scores .
Correlations between physiological parameters and clinical UPDRS sub-scores were calculated as described in Sharott et al, 2014. The percentage of STN units that were significantly locked to the midline EEG (recorded in all but one patient, n = 11) and the mean resultant vector length were correlated with the combined rigidity (subscore 22) and akinesia/bradykinesia (subscores 23-26 and 31) scores OFF medication across patients. For the percentage phase locked units, if a unit was significantly locked (Rayleigh test, p < 0.05) to any of the sub-bands within the low (12 - 24) or high (25 - 35) beta range (e.g. 10-15 Hz within the low beta range), that unit was counted as significantly locked for the wider range. For the mean resultant vector length, the maximum vector length for each unit was calculated for each frequency band, and these were averaged across all the units recorded in that patient. A bivariate regression model was used to calculate the Pearson correlation coefficient (R2) and associated p-value to evaluate the strength, significance and the best linear fit of the correlation. All Pearson correlation values had a Cook’s distance of > 1 suggesting that there were no outlying points (Cook, 1977).
Results
Cortico-subthalamic recordings in PD patients .
Our aim was to define the dynamics of beta frequency synchronisation between STN neurons and different cortical areas in patients with PD. To this end, we recorded surface EEG from 3 positions on the scalp using bipolar derivations with a central electrode, resulting in one midline and two lateral channels (defined as ipsilateral and contralateral to the coincident STN recording, Fig. 1A). In addition, in 7 patients we were able to record ECoG signals through a silver-ball inserted into the burr hole, approximately over the prefrontal cortex (Fig. 1B). STN-LFPs could be recorded from either macro-tips or the central micro-tip in 11 patients (Fig. 1C). STN single and multiunit activities were recorded from multiple tungsten microelectrodes (Fig. 1D). Synchronisation between STN unit activities and beta oscillations in the LFP and cortical activities could often be observed by eye (Fig. 2). The mean (33.67 spikes/s) firing rate of SUAs (Supplemental Fig. 5A) was in good agreement with well-isolated units in parkinsonian animals (Mallet et al., 2008b) and patients (Sharott et al., 2014). As reported previously (Weinberger et al., 2006; Sharott et al., 2014), the power spectra of subthalamic neuron spike trains and STN LFPs displayed peaks in the sub-beta and beta frequency ranges (Supplemental Fig. 5B, C). The mean power spectra of the ECoG and EEG only had a prominent peak in the sub-beta, but not beta range (Supplemental Fig. 5D). ECoG and midline-EEG channels appeared to have higher beta and lower gamma power than the lateral channels, but this was not statistically significant (Kruskal-Wallis ANOVA, p < 0.05).
LFPs underestimate subthalamic synchronisation with lateral cortical areas .

Frequency dependence of phase locking reveals topographically specific cortico-subthalamic time delays.
Overall, comparison of STN LFP and unit coherence with the cortical signals suggested that LFPs may underestimate the coupling of lateral cortical channels. In addition, peaks in coherence between motor cortical oscillations and STN units are often sharp (Shimamoto et al., 2013), rendering spectral estimates for defining the synchronisation parameters of STN units suboptimal. Thus, we further analysed the coupling of STN units to cortical beta oscillations by examining the phase locking of spike trains to 5 Hz bands within the broader beta band for STN LFPs and each cortical signal. Up to 20 to 25% of STN units recorded were phase locked to beta oscillations the STN LFP using this method, depending on the significance level used (Rayleigh Test, p < 0.01/0.05; Supplemental Fig. 7A). In agreement with coherence analysis (Figure 3B), the proportion of STN units that were significantly phase-locked at beta frequencies varied relatively little between cortical channels (Rayleigh Test, p < 0.05; Supplemental Fig. 7B-E). Having established that pool of STN neurons locked to oscillations at different frequencies across the beta range, we used the mean phase angle of locking to these narrow beta bands to recalculate the time delays between STN units and LFP, ECoG and EEG signals. This technique allowed us to pool neurons that synchronised to narrow frequency bands and to derive time delays from this pooled data (Supplemental Figs. 1-4).
Firstly, we analysed the relationship between STN units and the STN LFPs, where there should be little or no delay between the two signals. As expected, phase of maximal firing of STN units in relation to the STN LFP was constant, on the peak of the oscillation, across the whole beta range. This was evident at the level of single subjects (Fig. 4A) and across the whole data set (Fig. 4B). The mean phase distribution across the frequency ranges was flat, suggesting negligible time delay between the STN units and LFPs (Fig 4Bii, and see Supp. Fig 1Dii for comparison). Recomputing this delay for randomly selected, phase locked LFP-unit pairs across the beta frequency range led to a distribution of delays with no clear peaks (Fig. 4Biii), suggesting that there was no consistent delay across the sample of LFP-unit pairs (note that zero delay will not result as a peak as the real and shuffled data both give flat distributions). In contrast, phase histograms for the ECoG and midline EEG channels were approximately antiphase between low and high beta frequencies (Fig. 5Ai-ii, Bi-ii). This was the result of a linear shift in the phase of maximal firing as the centre beta frequency increased (Fig. 5Aiii, Biii). For the ipsilateral EEG, the phase histograms low and high frequencies looked similar, with maximal firing on the descending phase of the oscillation (Fig. 5C). Maximal firing to the contralateral EEG was shifted by around a quarter cycle between low and high beta frequencies (Fig. 5Di-ii). Plotting the firing probability across all frequencies revealed that this was caused by the maximal firing phase moving linearly over two cycles across the entire beta range (Fig. 5Diii), suggesting a longer lag between the contralateral EEG and STN units than for the frontal and midline channels.

Cortico-subthalamic beta dynamics are maintained during movement desynchronisation.  
During movement, STN neurons often display clear modulation of their firing rate (Rodriguez-Oroz et al., 2001; Hanson et al., 2012), whereas STN LFPs become less synchronised with cortical beta oscillations (Lalo et al., 2008; Talakoub et al., 2016). Next we examined whether there was any relationship between these single neuron and network level phenomena. To this end, patients were asked to perform periods of continuous extensions and flexion of the wrist contralateral to the recorded STN, which were monitored using EMG electrodes (Fig. 7A). Rate modulation of STN units could be readily observed online and individual units could both decrease and increase their firing rate during periods of continuous movement (Fig. 7B). A concurrent reduction in the phase locking of individual units could be observed during these during movement periods (Fig. 7C).
In order to quantify the relationship between the firing rate response and synchronisation to cortex, we calculated phase locking using a z-score that compared the real vector lengths for each unit to that of ISI shuffled surrogates. This method ensured that synchronisation could not be spuriously related to the number of spikes included in the analysis. For each unit we plotted this phase locking z-score at the frequency with maximum phase locking at rest, against the z-score at the same frequency during movement (Number of units. ECoG n = 25, Midline, n = 38, Ipsi. n = 41, Contra. n = 39. Fig. 8Ai-iv). Across all cortical channels, the phase locking of STN units was significantly lower during movement periods (P < 0.00000001, F = 107, Df = 1, effect of rest vs. movement, ANOVA). In addition, the phase locking of STN units varied significantly by cortical area (P = 0.003, F = 4.85, Df = 3, effect of rest vs. movement, ANOVA). Notably, the phase locking magnitude of STN units to the ECoG during movement was not significantly different to that of any of the cortical channels at rest and was significantly higher than that of the lateral channels (post-hoc Dunn-Sidak test). While the majority of STN units were less synchronised with the ECoG and midline EEG during the movement periods, around half of the units were still significantly phase locked (z-score > 2, ECoG = 60%, Midline = 46%. Fig. 8Aii-Bii). In contrast, relatively few STN units were significantly phase-locked to the ipsilateral and contralateral EEG during movement periods (z-score < 2; ipsilateral, 17%; contralateral, 21%), despite many units having high phase-locking values at rest (Fig. 8Aiii, Aiv). In line with these observations, the proportion of STN units significantly phase locked to frontal channels during movement was significantly higher than for lateral channels (Chi Square, p = 0.009).
Several studies have reported beta desynchronization in cortical signals during movement in PD patients (Litvak et al., 2012; Heinrichs-Graham et al., 2014; de Hemptinne et al., 2015; Rowland et al., 2015). If the amplitude of the beta signal was lower during movement periods, the drop in phase-locking strength of a given STN unit could be a consequence of a decrease in the signal-to-noise ratio of the cortical oscillation. However, we found that none of the cortical channels had a significant difference between rest and movement periods in the amplitude of the beta-frequency to which the STN units were phase-locked (Wilcoxon- sign rank test, p > 0.05). Moreover, there were no significant correlations between the percentage change in beta amplitude (at which the unit was phase-locked) and the percentage change or difference in z-scored vector length between rest and movement (Pearson Correlation, p > 0.05). While these analyses do not preclude cortical desynchronisation per se (which might be detected by triggering to movement onset, combining multiple frequency bands, over shorter timescales etc.), these analyses suggest that it did not underlie the changes in phase locking shown in Figure 8A.