Abstract
Oscillations in neural activity play a critical role in neural computation and communication. There is intriguing new evidence that the nonsinusoidal features of the oscillatory waveforms may inform underlying physiological and pathophysiological characteristics. Time-domain waveform analysis approaches stand in contrast to traditional Fourier-based methods, which alter or destroy subtle waveform features. Recently, it has been shown that the waveform features of oscillatory beta (13–30 Hz) events, a prominent motor cortical oscillation, may reflect near-synchronous excitatory synaptic inputs onto cortical pyramidal neurons. Here we analyze data from invasive human primary motor cortex (M1) recordings from patients with Parkinson's disease (PD) implanted with a deep brain stimulator (DBS) to test the hypothesis that the beta waveform becomes less sharp with DBS, suggesting that M1 input synchrony may be decreased. We find that, in PD, M1 beta oscillations have sharp, asymmetric, nonsinusoidal features, specifically asymmetries in the ratio between the sharpness of the beta peaks compared with the troughs. This waveform feature is nearly perfectly correlated with beta-high gamma phase-amplitude coupling (r = 0.94), a neural index previously shown to track PD-related motor deficit. Our results suggest that the pathophysiological beta generator is altered by DBS, smoothing out the beta waveform. This has implications not only for the interpretation of the physiological mechanism by which DBS reduces PD-related motor symptoms, but more broadly for our analytic toolkit in general. That is, the often-overlooked time-domain features of oscillatory waveforms may carry critical physiological information about neural processes and dynamics.
SIGNIFICANCE STATEMENT To better understand the neural basis of cognition and disease, we need to understand how groups of neurons interact to communicate with one another. For example, there is evidence that parkinsonian bradykinesia and rigidity may arise from an oversynchronization of afferents to the motor cortex, and that these symptoms are treatable using deep brain stimulation. Here we show that the waveform shape of beta (13–30 Hz) oscillations, which may reflect input synchrony onto the cortex, is altered by deep brain stimulation. This suggests that mechanistic inferences regarding physiological and pathophysiological neural communication may be made from the temporal dynamics of oscillatory waveform shape.
Introduction
Oscillatory activity is abundant in electrophysiological recordings in a variety of species, brain regions, and spatial scales (Engel et al., 2001; Buzsáki and Draguhn, 2004; Schnitzler and Gross, 2005; Buzsáki, 2006). These neural oscillations are critical for efficient communication within and across brain structures (Fries, 2005, 2015; Voytek and Knight, 2015). The vast majority of research examining the relationship between neural oscillations and cognition or disease uses Fourier-based methods to filter the signals into semiarbitrary frequency bands of interest. Although powerful, these methods inherently smooth the signal, altering or removing finer temporal features. This is especially striking given the mounting evidence that features of the oscillatory waveform, such as its general shape, sharpness, rise-to-decay symmetry, etc., may inform underlying physiology (Sherman et al., 2016; Cole and Voytek, 2017). In particular, some cortical oscillations in the beta (13–30 Hz) range have been shown to have characteristic waveforms that are shaper and steeper than a canonical sinusoid (Sherman et al., 2016). Computational modeling suggests these nonsinusoidal features may reflect the temporal synchrony of excitatory input currents onto cortical pyramidal neurons (Sherman et al., 2016), such that a sharper waveform reflects greater synchrony in synaptic input. That is, careful consideration of oscillatory waveform shape may provide critical clues about underlying neurophysiology and neural communication.
We examine the hypothesis that deep brain stimulation (DBS) treatment in patients with Parkinson's disease (PD) changes the waveform shape of beta oscillations in the primary motor cortex (M1), indicating that input synchrony to M1 may be altered. Previous reports have associated untreated PD with excessive beta synchronization in the basal ganglia (Brown, 2003). Additionally, recent evidence in M1 has shown increased phase-amplitude coupling (PAC) between beta phase and high gamma (50–200 Hz) amplitude in PD patients (de Hemptinne et al., 2013). This pathologically strong PAC is reduced during DBS (de Hemptinne et al., 2015). This is particularly intriguing, as high gamma amplitude has been shown to index local population firing rate (Mukamel et al., 2005; Ray et al., 2008; Manning et al., 2009; Miller et al., 2009a, 2014). Such exaggerated PAC in PD has been broadly interpreted to represent an oversynchronization of local spiking in the cortex (Voytek and Knight, 2015; de Hemptinne et al., 2015).
Here we observe that PD patients show a change in beta waveform shape during DBS, in line with the hypothesis that PD is associated with oversynchronized cortical inputs that are ameliorated via DBS. Additionally, the sharpness of a patient's beta oscillations predicted motor rigidity such that sharper beta (more synchronized inputs) was associated with greater rigidity. We extend the previous report showing decreased beta-high gamma PAC in PD patients during DBS, finding that beta waveform sharpness strongly tracks PAC across subjects. Although previous reports have shown that sharp transients can lead to “spurious” PAC (Kramer et al., 2008; Tort et al., 2013; He, 2014; Aru et al., 2015; Gerber et al., 2016; Lozano-Soldevilla et al., 2016), we emphasize here that these sharp beta features are not artifactual, in the sense that they are not the apparent result of mechanical or electrical noise. Rather, these features may provide novel insights into the physiological processes that generate them, with PAC being one method of detecting sharp waveforms that may reflect oversynchronized inputs. Because qualitatively different waveforms can underlie increased PAC (Vaz et al., 2017), we contend that a more thorough analysis of PAC should include time-domain characteristics in conjunction with spectral features.
Materials and Methods
Recordings analyzed in this report are the same as those analyzed in a previous report of 23 PD patients before and during DBS (de Hemptinne et al., 2015).
Data collection.
M1 recordings were obtained from 23 PD patients (20 male, 3 female) as previously described (de Hemptinne et al., 2013, 2015; Panov et al., 2017). Patients were recruited at the University of California at San Francisco or the San Francisco Veteran Affairs Medical Center. Patients were diagnosed with idiopathic PD with mild to moderate bradykinesia/rigidity and UPDRS III (Unified Parkinson's Disease Rating Scale part III) scores between 30 and 60. Patients underwent DBS implantation in the awake state and provided written informed consent. Patients were excluded if they had prominent tremor or peak-to-peak M1 LFP amplitude <50 μV. Data collection was approved by the institutional ethics committees and was in agreement with the Declaration of Helsinki.
DBS was as previously described (de Hemptinne et al., 2015). T2-weighted MRI was used to target the STN, with adjustments made based on movement-related spiking activity. Intraoperative computed tomography scans coregistered with preoperative MRI were used to confirm electrode placement of the DBS lead (model 3389 in 17 patients and 3387 in 6 patients; Medtronic). An analog neurostimulator (Medtronic model 3625) was used to set therapeutic stimulation parameters. Because optimal stimulation settings were not found before recording, an increased voltage (4 V) was used for stimulation between the motor territory of STN and its dorsal border. More details on patients and stimulation can be found in the previous reports (de Hemptinne et al., 2013, 2015).
A six-contact subdural electrocorticography strip was placed on the cortical surface using the burr hole for DBS lead placement. The target was the arm area of M1, 3 cm from midline and medial to the hand knob. Electrode contacts were platinum with a 4 mm total diameter, 2.3 mm exposed diameter, and 1 cm spacing between contacts (Ad-Tech). Correct placement of electrodes was confirmed using intraoperative computed tomography merged with preoperative MRI or lateral fluoroscopy. Additionally, physiological confirmation was obtained using median nerve stimulation (frequency = 2 Hz, pulse width = 200 μs, pulse train length = 160, amplitude = 25–40 mAmp) to evoke somatosensory potentials. The most posterior contact showing a negative N20 waveform was defined as the closest electrode to M1.
Antiparkinsonian medication was stopped 12 h before surgery. Data were collected 5–60 min after lead insertion to minimize the confounding effect of a temporary “microlesion” associated with lead insertion. In bilateral DBS implantation surgeries, brain activity was recorded on the second side implanted to allow more time between the cessation of propofol sedation and the start of electrocorticography recording. First, data were collected and evaluated after lead insertion, before any stimulation (“before DBS”). Second, data were collected when DBS was turned on for the first time (“during DBS”), before searching for optimal contact and stimulation parameters. Third, “after DBS” data were collected after DBS turned off for several minutes.
Recordings were collected using Microguide Pro (Alpha Omega) or the Guideline 4000 customized clinical recording system (FHC) with a sampling rate between 1000 and 3000 Hz. A needle electrode in the scalp was used as the ground. Signals were bandpass filtered 1–500 Hz and amplified 7000×. While the analyzed data were collected, subjects were relaxing with eyes open, fixating on a point ∼1 m away. The first 30 s of data without obvious electrical noise or movement were selected for analyses.
Data preprocessing.
Data were processed in the same way as previously reported (de Hemptinne et al., 2013, 2015). Recordings were referenced using a bipolar montage in which each channel was referenced to the immediately anterior channel. Data were downsampled to 1000 Hz. Line noise was removed with a notch filter between 58 Hz and 62 Hz (Butterworth, order = 3). Additionally, in this paper, all signals were cleaned of high-frequency artifacts by low-pass filtering at 200 Hz (FIR filter, window method, order = 250 ms) and applying narrowband notch filters (Butterworth, order = 3) to remove any sharp peaks in power spectra >80 Hz caused by DBS stimulation and electronic noise. The same filters were applied to all recordings from the same subject (before DBS and during DBS).
Oscillation sharpness.
The procedure for estimating the sharpness of beta oscillations is illustrated in Figure 1A, B and begins with finding oscillatory peaks and troughs. First, the raw voltage trace is bandpass filtered using an FIR filter (window method, cutoff frequencies = 13 and 30 Hz, order = 231 ms). Time points of rising and falling zero-crossings are identified. Returning to the raw signal, the time point of maximal voltage between a rising zero-crossing and a subsequent falling zero-crossing is defined as the peak. Similarly, the time point of minimal voltage between a falling zero-crossing and a subsequent rising zero-crossing is defined as a trough.
Sharpness of each peak and trough is defined in Equation 1 in which Vpeak is the voltage at the oscillatory peak and Vpeak−5 ms and Vpeak+5 ms correspond to the voltage 5 ms before and after the peak, respectively, as follows:
Trough sharpness is calculated in the same manner. Intuitively, extrema sharpness increases as the absolute voltage difference between the extrema and the surrounding time points increases. Although 5 ms is chosen as the temporal width to quantify sharpness, results are similar for other choices of sharpness width (data not shown).
The mean sharpness throughout a recording is calculated for peaks and troughs separately (∼600 each in each 30 s recording). The reported sharpness ratio metric is the ratio of the sharpness of the two extrema, such that the ratio is always >1 as in Equation 2 as follows:
This ratio accounts for differences in the amplitudes across subjects that can be caused by variance in electrode conductance. In Figure 2, this ratio is not fixed to be >1, but rather mean peak sharpness is divided by mean trough sharpness.
Oscillation steepness.
Steepness of the rises and decays of each oscillation is calculated using the peaks and troughs identified as described above. For rise steepness, the time series of interest is the raw voltage between each trough and the subsequent peak as in Equation 3 as follows:
The maximum value of the instantaneous first derivative (numpy.diff) of this signal is defined as the steepness for this single rise period. The steepness of decay periods is similarly calculated using the raw voltage time series between each peak and the subsequent trough. The average rise and decay steepness are calculated for each recording across all rise and decay periods. A steepness ratio is then calculated for each recording by dividing the average steepness for all rises/decays by the average steepness for all decays/rises (such that the ratio is >1) as in Equation 4 (below). In Figure 2, this ratio is not fixed to be >1, but rather rise steepness is divided by decay steepness.
Specific code for the novel metrics introduced, along with tutorials, is available at https://github.com/voytekresearch/misshapen.
Spectral analysis.
The power spectrum of a signal is calculated by squaring the absolute value of its Fourier transform (numpy.fft.fft). Beta power (see Fig. 5) is calculated for each 30 s signal by summing the coefficients in the power spectrum in the beta frequency range (13–30 Hz). Spectrograms of signals (see Fig. 4E) are calculated using a continuous wavelet transform (7 cycle complex Morlet wavelets 1–200 Hz with a 1 Hz step size, scipy.signal.morlet with s = 0.5). Beta frequency for each recording is estimated by counting the number of peaks identified as described above and dividing by the duration of the recording.
The beta and high gamma components of the signal are obtained by using FIR filters (window method, scipy.signal.firwin). Two-pass zero-phase filtering is performed with these filters (scipy.signal.filtfilt). Filter orders are chosen to obtain both a desirable frequency response and reasonable temporal resolution. Beta bandpass filters have cutoff frequencies of 13 and 30 Hz and an order of 231 ms. High gamma bandpass filters have cutoff frequencies of 50 and 200 Hz and an order of 240 ms. Amplitude of the high gamma component (see Fig. 6A,B) is calculated by the magnitude of its Hilbert transform. The high gamma amplitude is sampled at the peaks and troughs in each beta oscillation (see Fig. 4B).
Statistical PAC is estimated using the normalized modulation index metric (Özkurt and Schnitzler, 2011). PAC was estimated using the open-source package, pacpy, available at https://github.com/voytekresearch/pacpy. The high gamma amplitude is calculated as described above. The beta phase is similarly calculated by the angle of the Hilbert transform of this component. Results are similar when applying four alternative metrics of statistical PAC (Canolty et al., 2006; Penny et al., 2008; Tort et al., 2010) or when a comodulogram method is used as in the previous report (de Hemptinne et al., 2015) (data not shown). The preferred phase of coupling for each recording (see Fig. 3G) is determined by the angle of the circular sum of the instantaneous beta phase, weighted by the instantaneous high gamma amplitude.
Comodulograms (see Figs. 3A,B, 6C,D) are calculated similar to the previous report (de Hemptinne et al., 2015). The frequency range for the phase-providing oscillation ranges from 6 to 40 Hz with a step size and bandwidth of 2 Hz. For the amplitude-providing oscillation, the frequency ranges from 20 to 200 Hz with a step size and bandwidth of 4 Hz. At each frequency step, a 7 cycle complex Morlet wavelet is used as a bandpass filter, and either the instantaneous angle or magnitude of the filtered signal is calculated to extract phase or amplitude time series, respectively. The modulation index method (Tort et al., 2010) is used to quantify statistical PAC between each combination of filtered signals in the comodulogram.
Statistical PAC histograms (see Fig. 3D–F) are calculated from the beta phase time series and high gamma amplitude time series used to calculate statistical PAC. The high gamma amplitude is then averaged across 10 equally sized bins of beta phase. A similar phase estimate as described previously (Siapas et al., 2005; Belluscio et al., 2012; Trimper et al., 2014) is applied to account for the nonsinusoidal shape of neural oscillations. Time points of peaks and troughs are identified as in Oscillation sharpness. The time points of these extrema are then used to linearly interpolate a theoretical phase value for each sample.
Analysis of individual beta cycles is performed by separating the cycles trough-to-trough or peak-to-peak using the extrema found as described above. The sharpness of each cycle is measured, and they are split into five groups based on this value to compare extrema-locked high gamma amplitude (see Fig. 4C) or statistical PAC (see Fig. 4D). In Figure 5C, the sharpness ratio between each trough and subsequent peak is correlated with the voltage difference between these two extrema.
To visualize potential coupling between beta and gamma oscillations, event-related averages of the raw data (see Fig. 4F) were calculated by triggering on the peaks of the high gamma component (obtained with the filter described above). High gamma peaks were only used as triggers if the amplitude at that time point was in the top 10th percentile, to limit triggering during periods with negligible high gamma power in which the peak phase is essentially random.
Generation of canonical asynchronous PAC.
A canonical signal with beta-high gamma asynchronous PAC (see Fig. 6B) is simulated in three steps: (1) a beta oscillation is simulated by bandpass filtering 30 s of white noise with a bandpass FIR filter (13–30 Hz); (2) high gamma is similarly simulated by bandpass filtering 30 s of white noise with a 50–200 Hz bandpass FIR filter (the amplitude of the high gamma is then modulated by the beta phase by multiplying its time series by 1 [− abs(φβ)/π] and scaling by 0.03 to decrease its amplitude relative to the simulated beta oscillation); and (3) the beta oscillation and high gamma components are added together.
Statistics.
The Scipy package (version 0.16.0) in Python (version 2.7) is used for all statistical analysis. Unless indicated otherwise, all correlations are Pearson and all p values are two-tailed.
Results
M1 beta oscillations in PD are flattened by DBS
To quantify beta waveform shape, peaks and troughs were first identified throughout each 30 s recording (Fig. 1A; see Materials and Methods). After locating the extrema, the sharpness of each peak and trough was quantified (Fig. 1B; see Materials and Methods). Importantly, the metric we use for extrema sharpness is higher for oscillations with sharp transients (e.g., sawtooth waves) compared with sinusoids. The symmetry in the sharpness of oscillatory peaks and troughs is visualized by the degree of overlap in the distributions of sharpness over all extrema in the signal (Fig. 1C,D) and is quantified as the “sharpness ratio.” We find that DBS treatment decreases sharpness ratio in PD patients (Fig. 1E; paired t test, t(22) = 2.5, p = 0.019), and thus, asymmetrically affects peak and trough sharpness. Furthermore, PD patients' clinical rigidity scores before DBS positively correlate with sharpness ratio (Fig. 1F; Spearman r = 0.54; n = 23; p = 0.014). However, the changes in rigidity score and sharpness ratio do not correlate with each other (Spearman r = 0.17; n = 23; p = 0.48).
Waveform shape of M1 beta oscillation changes in PD. A, Estimation of beta oscillatory extrema. Row 1, a raw voltage signal; Row 2, the raw voltage is bandpass-filtered in the beta frequency range to identify rising phases (yellow dots) and falling phases (green dots); Row 3, in between each rise and decay, peaks (purple circles) and troughs (teal circles) are identified in the raw voltage signal. B, Sharpness of each extrema is calculated as shown by the formulas on the right. The extrema sharpness is the average difference between the voltage at the extrema and the voltage at 5 ms before and after the extrema. These time points are indicated by teal vertical lines around the trough and purple vertical lines around the peak. The voltage ±5 ms around the trough is around the same as the trough voltage (blue triangles and circle), and the voltage ±5 ms around the peak changes relatively more (purple triangles, relative to purple circle). C, D, Distributions of peak and trough sharpness across a 30 s recording for a PD patient (C) before DBS and (D) during DBS. E, The sharpness ratio (between peak sharpness and trough sharpness) is decreased by DBS. F, Clinical rigidity scores are positively correlated with sharpness ratio in PD patients before DBS. E, F, Each dot indicates 1 patient. *p < 0.05.
Differences in sharpness ratio can be caused by increases and/or decreases in the sharpness of one or both extrema. Here, we find that the average sharpness of peaks and troughs decreases with DBS application (paired t test, t(45) = 2.4, p = 0.027). Additionally, the extrema sharpness itself does not correlate with the clinical rigidity score (Spearman correlation; r = −0.17; n = 23; p = 0.49), indicating that it is the ratio of peak to trough sharpness that is changed after DBS, and not their sharpness in general.
M1 beta has a consistent sawtooth-like waveform
To characterize the waveform of M1 beta, we additionally compute a steepness ratio, similar to the sharpness ratio, across subjects. The steepness ratio quantifies the asymmetry between the steepness of the rise and the steepness of the decay period within a beta cycle (see Materials and Methods). Across recordings, steepness ratio is strongly correlated with sharpness ratio (Fig. 2A; r = 0.84, p < 10−12). That is, beta oscillations whose peaks are sharper than their troughs have rise phases that are steeper than their decay phases (Fig. 2B, top right). In contrast, oscillations whose troughs are sharper than their peaks have decay phases that are steeper than their rise phases (Fig. 2B, bottom left). Example recordings from 2 PD patients before DBS show these sawtooth-like waveform shapes in raw data (Fig. 2C). The neural dynamics in M1 seem to produce field potentials that have a consistent sawtooth-like waveform in which the extremum following a steep voltage change is sharper than the extremum preceding the steep voltage change.
M1 beta oscillations have a consistent sawtooth shape. A, Positive correlation between the relative sharpness beta oscillatory peaks and the relative steepness of beta voltage rises. Each marker represents the average sharpness and steepness in a 30 s recording from a PD patient either before DBS (blue) or during DBS (red). B, Schematic voltage traces corresponding to each quadrant of A. M1 beta falls in the dark gray quadrants (quadrants I and III) of this two-dimensional space. C, Example voltage traces from 2 patients corresponding to the sawtooth waveforms in quadrants I and III. These waveform shapes are consistent with the dark gray, but not the light gray, sawtooth shapes in B.
Beta oscillation sharpness increases statistical PAC
Previous work showed that beta-high gamma PAC is also changed by DBS application (de Hemptinne et al., 2015), and additionally showed that PAC in PD patients was higher relative to epilepsy and cervical dystonia patients (de Hemptinne et al., 2013). In the following, we extend these findings by relating PAC to sharpness ratio. This is interesting because the relation of these past results to waveform shape can lead to novel hypotheses regarding the relevant pathophysiological processes in PD (Cole and Voytek, 2017). Because of recent reports (Gerber et al., 2016; Lozano-Soldevilla et al., 2016; Vaz et al., 2017), we make it explicit that our measure of PAC is agnostic to the type of processes that generates it; we will refer to it as statistical PAC in the remainder. Importantly, we reiterate that this by no means indicates that the PAC we quantify is artifactual (i.e., due to nonbrain sources), only that it can reflect multiple generative mechanisms within the brain. The statistical relationship between these two frequency bands is true, although the common physiological interpretation of two coupled processes differs from the interpretation as changes in waveform sharpness.
First, we show that M1 beta-high gamma statistical PAC was decreased during DBS application, following previous reports. The comodulograms for 1 example patient (Fig. 3A,B), shows higher coupling between beta phase and broadband high gamma amplitude before DBS compared with during DBS. In general, DBS application to PD patients decreases the statistical PAC (Fig. 3C; paired t test, t(45) = 2.7, p = 0.013), as previously reported. We newly observe that, in general, high gamma amplitude is specifically coupled to the peaks and troughs of the beta oscillations, as opposed to nonextrema phases of the beta cycle (such as the rise or decay periods of the oscillation). This is shown for 3 example PD patients before DBS, with the strongest high gamma at the peak (Fig. 3D), trough (Fig. 3E), or both (Fig. 3F), and for all PD patients in Figure 3G. This pattern of increased PAC at peaks and troughs is statistically observed as a positive correlation between the modulation index and the absolute value of the cosine of the preferred phase (Spearman correlation, r = 0.58, n = 23, p = 0.003).
Estimates of phase-amplitude coupling in PD. A, Comodulogram showing estimated beta-high gamma PAC in an example PD patient before DBS. B, The estimated PAC decreases with DBS application. C, Estimated PAC decreases with DBS application. D–F, Distributions of high gamma amplitude as a function of beta phase for 3 example PD patients before DBS showing increased high gamma amplitude at (D) the peak of the beta oscillation, (E) trough of the beta oscillation, and (F) both extrema of the beta oscillation. G, High estimated PAC is observed specifically in recordings in which the high gamma amplitude is increased around the peak (phase = 0) or trough (phase = −π, π) C, G, Each dot indicates 1 patient. *p < 0.05.
We now extend the above by presenting a converging sequence of results that shows the extrema sharpness of beta waveforms is strongly related to statistical PAC (Fig. 4), providing strong evidence that sharpness is the main contributor to phase-locked high gamma. First, we find that the sharpness ratio of beta oscillations is strongly correlated with beta-high gamma statistical PAC across PD patients before DBS (Fig. 4A; r = 0.94, p < 10−10) and during DBS (r = 0.89, p < 10−7). The remaining results are illustrated for one representative subject before DBS (Fig. 4B–F), and similar results were obtained in all other recordings. Second, we find a positive correlation between extrema sharpness and high gamma amplitude at that extrema (Fig. 4B). The correlation between these features is similar for both peaks and troughs. Third, providing a time-resolved view of the above, the strength of high gamma amplitude increases as a function of extrema sharpness (Fig. 4C). Fourth, the strength of phase-locking of high gamma amplitude across the beta cycle increases as a function of extrema sharpness (Fig. 4D). This is known to be the case when PAC is driven by sharp temporal features (He, 2014). The asymmetry in high gamma amplitude around the beta peak (Phase 0) indicates stronger high gamma amplitude in the rise phase of the beta oscillation relative to its decay phase. This follows from the sawtooth-like shape of the beta oscillations (Fig. 2B, quadrants 1 and 3), in that high-frequency sinusoids are necessary to reconstruct steep voltage rises. Fifth, we show a time-frequency representation of amplitude (Fig. 4E) for an example burst of beta oscillations (white trace, peaks numbered) of varying sharpness (values printed above plot). We observe that the high-frequency amplitude is strongest (highest PAC) for the sharper peaks (peaks 2, 3, 4, and 6), and weakest (lowest PAC) at the smoother peaks (peaks 1 and 5). Sixth, and finally, we show an event-related average triggered on the peaks of the high gamma component (see Materials and Methods). This is a common type of visualization for PAC. Typically, in the case of two coupled oscillations, this reveals a waveform in which the high-frequency oscillation appears ∼t = 0, together with the low-frequency oscillation at an offset matching the phase of the statistical PAC (Tort et al., 2013). However, in Figure 4F, we do not observe a high-frequency oscillation, but rather see a sharp voltage transient, as would be expected from a sawtooth-like waveform (Tort et al., 2013).
Sharp beta oscillations yield phase-locked high gamma amplitude. A, Sharpness ratio is highly correlated with estimated PAC across all subjects before DBS. B–F, Illustration of the relationship between extrema sharpness and phase-locked high gamma amplitude for an example PD patient before DBS. B, Extrema sharpness is positively correlated with the high gamma amplitude at the extrema. C, High gamma amplitude is locked to the peak of beta oscillations specifically for the sharpest beta oscillations. Beta cycles were split into five groups based on the peak sharpness. Line color varies from the least sharp beta cycles (light gray) to the sharpest beta cycles (black). D, The distribution of high gamma amplitude as a function of beta phase is uniform for beta cycles with low sharpness (light gray) and nonuniform for sharp beta cycles. E, Spectrogram of a 400 ms period of raw voltage recording exhibiting beta oscillations. The sharpest beta oscillations (peaks 2, 3, 4, and 6) have increased amplitude in the high gamma frequency range compared with the cycles with lower sharpness (peaks 1 and 5). F, Event-related average of raw voltage triggered at the peaks of high gamma oscillations (with the top 10% of amplitude).
Caveats of the relation between waveform sharpness and statistical PAC
A positive correlation between sharpness ratio and statistical PAC could be caused by sharp, phase-locked high gamma rhythms, instead of the sharpness of beta oscillations highlighted above. To gain more certainty that sharpness came from beta instead of high gamma, sharpness ratio was recalculated on data that was low-pass filtered at 50 Hz (FIR filter, window method, order = 250 ms). This removes the high gamma component used to estimate PAC (Fig. 5A). This new sharpness ratio measure (after low-pass filtering) is still correlated with the statistical PAC computed before low-pass filtering (Fig. 5B; r = 0.80, p < 10−5). Some decrease in correlation strength (in this case, r = 0.94 to 0.80) is expected because sharp extrema are (by definition) broadband in their spectrum and are necessarily flattened by the 50 Hz low-pass filter. Crucially, although the beta waveform is only mildly affected by the low-pass filter, the high gamma fluctuations are suppressed to the level of noise floor of the amplifier. This makes it extremely unlikely that any residual high gamma fluctuations cause a correlation, let alone one of r = 0.80. As such, the correlation provides strong evidence that the sharpness measures are not influenced by high gamma fluctuations, but instead originated from the shape of the beta waveform.
Power spectral properties of motor cortical recordings and relation to estimated phase-amplitude coupling. A, A 50 Hz lowpass filter was applied to attenuate high gamma power. B, When extrema sharpness was recalculated on these lowpass-filtered signals, there remained a strong positive correlation with statistical PAC in the original signal. C, Histogram of correlation coefficients between sharpness ratio and amplitude across all beta cycles in a recording. Each count is one 30 s recording. D, There is a positive correlation between estimated PAC and beta power. E, A strong positive correlation between sharpness ratio and estimated PAC (modulation index) remains even after beta power is regressed out. F, Average power spectra across all subjects before DBS (black) and during DBS (blue). Shaded region represents SEM. B, D, E, Each dot indicates one 30 s recording.
In our approach, the sharpness is calculated at each beta cycle throughout the recording, even if no beta oscillation is evident in the raw voltage. One concern with this approach is that noise in a regimen of low beta amplitude could cause spurious results in waveform analysis. To test for this potential confound, we restrict our analysis to only the cycles in the top 10% of beta amplitude for each recording. The main results still hold under this control analysis: (1) there is a decrease in sharpness ratio with DBS application (paired t test, t(22) = 2.6, p = 0.015); and (2) there is a strong correlation between sharpness ratio and PAC (r = 0.84, p < 10−6).
Because sharpness ratio is calculated using the raw electrophysiological signal of each extremum and two samples (sampled at 1 kHz) around it, it is worth noting that 90% of the variance in statistical PAC is captured by only using ∼12% of the data (6 samples per period, which contains ∼50 samples). However, sharpness ratio is not the only dimension of shape that correlates with statistical PAC, as many other waveform features determine the precise Fourier decomposition (such as steepness ratio). Sharpness ratio and steepness ratio both individually correlate with statistical PAC after holding the other metric constant (partial correlations, sharpness ratio: r = 0.70, p = 0.0002; steepness ratio: r = 0.44, p = 0.034) and together explain 93% of the variance in statistical PAC of PD patients before DBS.
Another potential caveat in the analysis of waveform sharpness is that increases in beta frequency could underlie increases in the sharpness, as a shorter cycle of the same amplitude is necessarily sharper. However, this is likely not the case, as we find no effect of DBS on the beta frequency in PD patients (paired t test, t(22) = −0.84, p = 0.41; see Materials and Methods). Power does not confound the sharpness ratio measure because multiplicative scaling of a waveform to increase its power would have no effect on the sharpness ratio.
Beta oscillation asymmetry correlates with power
During visual inspection of the raw data, high-power beta oscillations seemed to be more asymmetric (i.e., higher sharpness ratio) compared with low-power beta oscillations. If true, this trend could provide additional insight into the physiological process generating beta oscillations. We quantified this trend by correlating the sharpness ratio of an individual cycle with the peak-to-trough amplitude of that cycle. To avoid confounding the results with periods of no beta oscillatory activity, only cycles in the top 10th percentile of peak-to-trough amplitude for each recording were analyzed. Across recordings before and during DBS, there is a consistent subject-by-subject positive correlation between the amplitude of a cycle and its sharpness ratio (Fig. 5C; one-sample t test, t(45) = 6.10, p < 10−6).
Because beta power is positively correlated with both sharpness ratio (Fig. 5C) and statistical PAC (Fig. 5D; r = 0.37, p = 0.011), it could be the case that the correlation between waveform shape and statistical PAC is merely due to a common power bias. However, this is likely not the case, as sharpness ratio is still correlated with statistical PAC after holding beta power constant (Fig. 5E; partial correlation, r = 0.87, p < 10−7). Additionally, we note that the average power spectra are similar before and during DBS (Fig. 5F), and so it is unlikely that the changes in waveform shape and statistical PAC with DBS are merely due to changes in narrowband power.
Discussion
Increased sharpness ratio in motor cortical beta oscillations in Parkinson's Disease
Although beta oscillations are a normal feature of the basal ganglia-thalamocortical loop, PD is associated with excessive neuronal synchronization in the beta frequency band (Kühn et al., 2005; Moran et al., 2008). Despite an established relationship between beta band neuronal synchronization and PD, the physiological mechanism causing PD-related motor dysfunction has been unclear. However, the recent reports showing that statistical PAC is pathologically strong in PD hypothesize that this overcoupling impairs information flow and in turn causes motor dysfunction (de Hemptinne et al., 2013). Specifically, PD is associated with an increase in beta-high gamma statistical PAC in M1 (de Hemptinne et al., 2013, 2015; Kondylis et al., 2016). Similarly, elevated alpha-high gamma statistical PAC has been observed in the sensorimotor cortex of essential tremor patients (Kondylis et al., 2016). These observations have been critical given that oscillatory coupling, and specifically PAC, is thought to be crucial for the control of information flow in human cortex (Fries, 2005; Canolty et al., 2006; Peterson and Voytek, 2015; Voytek and Knight, 2015; Voytek et al., 2015).
We confirm and extend the previous report of exaggerated statistical PAC in PD by analyzing oscillation shape, given the recent reports that oscillation waveform shape may provide information about their physiological generators (Sherman et al., 2016; Cole and Voytek, 2017). Specifically, recent modeling work suggests that cortical beta waveform sharpness may index the degree of synchrony of input onto cortical pyramidal cells (Sherman et al., 2016). Therefore, here we leveraged analytic approaches on time series to test the hypothesis that beta waveform sharpness is decreased by DBS, suggesting that the synchrony of synaptic input to M1 is decreased. For M1, synchronous bursts of distal excitatory input may reflect the previously reported increase in beta synchrony in the basal ganglia (Sharott et al., 2005; Mallet et al., 2008a, b; Devergnas et al., 2014). Just as exaggerated PAC has been hypothesized to overwhelm neural processing, too much synchrony consumes an excessive amount of “neural bandwidth” (Reyes, 2003; Rossant et al., 2011; Brette, 2012; Börgers et al., 2014), hampering, rather than enhancing, neural communications. Given how potent synchronous activity can be in driving downstream population spiking (Reyes, 2003; Wang et al., 2010; Rossant et al., 2011), we suggest that excessive synchrony in M1 may be key to understanding PD pathology. Therefore, our results support previous hypotheses that DBS reduces PD-related movement symptoms by decorrelating the excessively synchronized neural activity in the basal ganglia-thalamocortical loop (Rubin and Terman, 2004; Moran et al., 2012; Wilson, 2013; Voytek and Knight, 2015). Notably, the reported change in waveform shape with DBS leads to a similar interpretation as the decrease in estimated PAC. However, the critical difference between these two interpretations is that the former invokes a modification of a single process whereas the latter invokes a decrease in coupling between two processes.
We note that the higher-power beta oscillations tend to be the most asymmetric in terms of extrema sharpness (Fig. 5C). If the beta waveform reflects the strength and synchrony of the afferent synaptic currents, as reasoned above, this trend may reflect that beta periods with the strongest synaptic input also have the most synchronous input.
The current description of nonsinusoidal waveforms in motor cortex is reminiscent of previously described mu rhythms. Mu rhythms are arch-shaped (i.e., high sharpness ratio, unity steepness ratio) oscillations ∼10 Hz (Tiihonen et al., 1989; Pfurtscheller et al., 1997). Similarly, the present analysis observed more extreme sharpness ratios than steepness ratios (Fig. 2). Therefore, it may be tempting to relate the oscillations in the current study to the mu rhythm. However, the oscillations here are at a beta frequency, whereas mu rhythms characteristically oscillate at an alpha frequency. Therefore, the oscillations in the present study are likely related to previous reports of sensorimotor beta rhythms (e.g., Miller et al., 2012) and distinct from previous reports of mu rhythms.
Nonsinusoidal oscillations can underlie statistical PAC
Cross-frequency coupling (CFC) analysis techniques quantify interactions between neural oscillations. One example of CFC is PAC, in which the phase of an oscillation is correlated with the amplitude of a higher-frequency oscillation or broadband high-frequency activity (Canolty et al., 2006). The degree of interaction between these two frequency components has been estimated using PAC metrics, often with the assumption that low-frequency oscillatory phase organizes neuronal cell assembly spiking (Canolty et al., 2006; Lisman and Jensen, 2013; de Hemptinne et al., 2013, 2015; Voytek et al., 2015; Watrous et al., 2015). In this view, there are two separate interacting neural processes: a low-frequency oscillator that modulates local spiking probability.
However, it is helpful to distinguish several biophysical processes that can contribute to statistical PAC: coupling can occur between two oscillatory processes; alternatively, a low-frequency oscillator associated with synaptic currents (Miller et al., 2009b; Okun et al., 2010; Buzsáki et al., 2012; Einevoll et al., 2013; Mazzoni et al., 2015) may be coupled to high gamma activity, associated with asynchronous spiking activity (Ray et al., 2008; Manning et al., 2009; Miller et al., 2009a, 2014). Furthermore, statistical PAC is increased by nonsinusoidal waveforms, as previously been shown in synthetic data (Kramer et al., 2008; Lozano-Soldevilla et al., 2016), high-voltage spindles (Tort et al., 2013), and cortical alpha oscillations (Lozano-Soldevilla et al., 2016). These previous reports, as well as a recent review (Aru et al., 2015) have offered suggestions for determining the nature of PAC, which can be diverse across the human cortex, simultaneously exhibiting both oscillatory and nonsinusoidal modes (Vaz et al., 2017). We here emphasize the importance of visually inspecting the raw time series to judge whether statistical PAC is driven by sharp voltage changes, as shown here, by true oscillation-to-oscillation coupling, or by sustained asynchronous activity during particular low-frequency oscillatory phases.
The sharp waveforms reported here produce similar statistical PAC results compared with canonical asynchronous PAC waveforms (Fig. 6). The difference between the temporal dynamics in these waveforms (Fig. 6A,B) suggests a difference in the underlying neural activity, which is not differentiated by statistical PAC. While high gamma amplitude correlates with local population spiking activity (Ray et al., 2008; Manning et al., 2009; Miller et al., 2009a, 2014), the magnitude of previously reported high gamma changes are low, on the order of a few microvolts. In contrast, the apparent high gamma resulting from the sharp time-domain deflections seen here are an order of magnitude stronger, >100 μV in some cases (Fig. 6A). This extreme difference in magnitude suggests that the underlying phenomena of high gamma may not be the same in these two cases.
Sharp beta oscillations produce similar cross-frequency coupling estimates compared with canonical phase-amplitude coupling. A, Raw voltage recording in a PD patient before DBS (black) overlaid with the high gamma component (thick red line) and its analytic amplitude (thin red line). B, Same format as A with a simulated signal with canonical PAC between the phase of beta oscillations and high gamma amplitude. C, D, Comodulograms showing similar estimated PAC between beta phase and high gamma amplitude in both the (C) M1 recording and (D) canonical PAC signal.
In conclusion, we offer a new perspective in which sharp voltage transients within an oscillation may carry physiological information. These often-overlooked temporal domain features should complement spectral CFC analyses. When used in conjunction, temporal domain analysis can offer novel insights into the biophysical processes generating the statistical PAC. Here, we demonstrate that the sharpness of motor cortical beta waveforms, previously shown to reflect synchronous input, are decreased with DBS treatment of PD. Furthermore, the sharpness ratio measure positively correlates with clinical rigidity measures. This synchrony interpretation offers new insight into the pathophysiology of PD, serving as an example of how a combination of both spectral and temporal analyses may be useful in extracting critical information from electrophysiological signals.
Footnotes
S.R.C. was supported by the National Science Foundation Graduate Research Fellowship Program. B.V. was supported by the University of California–San Diego, Qualcomm Institute, California Institute for Telecommunications and Information Technology, Strategic Research Opportunities Program, and a Sloan Research Fellowship. B.V. and R.v.d.M. were further supported by National Institutes of Health Grant MH095984 to B.R. Postle. Data collection by P.A.S. and C.d.H. was supported by the Michael J. Fox Foundation and National Institutes of Health Grant R01 NS069779. We thank Tom Donoghue, Richard Gao, Torben Noto, Brad Postle, Tammy Tran, and Andrew Watrous for invaluable discussion and comments.
The authors declare no competing financial interests.
- Correspondence should be addressed to Scott R. Cole, University of California–San Diego, 9500 Gilman Drive, La Jolla, CA 92093-0515. scott.cole0{at}gmail.com
