Differential contributions of subthalamic beta rhythms and 1/f broadband activity to motor symptoms in Parkinson’s disease

Excessive beta oscillatory activity in the subthalamic nucleus (STN) is linked to Parkinson’s Disease (PD) motor symptoms. However, previous works have been inconsistent regarding the functional role of beta activity in untreated Parkinsonian states, questioning such role. We hypothesized that this inconsistency is due to the influence of electrophysiological broadband activity —a neurophysiological indicator of synaptic excitation/inhibition ratio— that could confound measurements of beta activity in STN recordings. Here we propose a data-driven, automatic and individualized mathematical model that disentangles beta activity and 1/f broadband activity in the STN power spectrum, and investigate the link between these individual components and motor symptoms in thirteen Parkinsonian patients. We show, using both modeled and actual data, how beta oscillatory activity significantly correlates with motor symptoms (bradykinesia and rigidity) only when broadband activity is not considered in the biomarker estimations, providing solid evidence that oscillatory beta activity does correlate with motor symptoms in untreated PD states as well as the significant impact of broadband activity. These findings emphasize the importance of data-driven models and the identification of better biomarkers for characterizing symptom severity and closed-loop applications.


Subjects and experimental protocol
LFPs were recorded between 36 and 48 hours after STN DBS bilateral implantation in thirteen subjects using a 16-channel amplifier (gUsbAmp, g.tec, Austria) and sampled at 512 Hz. All

Surgical procedure
Long-acting dopaminergic medication was withdrawn 48 hours prior and short-acting medication was withdrawn 12-20 hours prior to off-medication pre-operative testing and to the DBS lead implantation procedure, always performed by the same surgeon. The awake surgical procedure stereotactically implanted a 3389 Medtronic electrode in the motor part of the STN. Direct surgical planning was performed with the Medtronic stealth station, based on a CT scan obtained with the CRW stereotactic frame, and fused with two 3 Tesla MRI sequences (MPRAGE with Gadolinium and space T2), obtained prior to surgery. The precise trajectory of the implanted electrode was simulated and its exact coordinates were calculated in relation with the position of the stereotactic frame. During the surgery, a cannula was implanted towards the STN following the calculation of the trajectory and remained in place until the end of the surgery. Through this cannula, a microelectrode was lowered toward the target in the STN and the border of the motor STN was visualized as the transition from slow spiking neurons to fast spiking neurons. After a first phase of microrecording, macrostimulation was performed under neurological clinical assessment; improvement of rigidity was expected in case of correct placement. Adjustment of the electrode placement through a second track was eventually performed to optimize the clinical response to stimulation. The electrode location was checked during the surgery with the O-Arm (Medtronic).
An externalized extension cable was connected to the distal part of the electrode and tunneled posterior to the skin incision to record the brain activity for a period of 3 days. Internalization of the electrode and connection of the electrode to an internal stimulator were performed at day 4 under general anesthesia. Once the definitive electrode was placed through the same cannula to the definitive target, a second CT scan was obtained after having implanted the definitive electrode.
Finally, the electrode was fixed with cement on the skull.
In order to assess the efficacy of stimulation, we calculated the clinical score during OFF-versus ON-stimulation. Results showed that the clinical scores were significantly better during ONstimulation than during OFF-stimulation across patients and hemispheres (meanON= 2.6±2.2, meanOFF= 5.2±3.9; p=0.0001, signed rank test). The improvement in clinical score from OFF-to ON-stimulation was calculated as follows: The average improvement across patients and hemisphere was 43±7%, similar to those obtained by previous works 1 .

Power spectrum modeling
Six minutes of bipolar LFP activity were derived from two contacts, selected based on the electrodes chosen by the neurosurgeon for DBS stimulation. LFPs were first high-pass filtered at 1 Hz using a second order Butterworth filter. Artifact were removed using an automatic threshold approach (data with activity higher than 50µV in absolute value), followed by visual inspection and rejection of contaminated portions of data. Those windows with abnormal activations (higher than 50 µV for the automatic approach; or marked as noisy by visual inspection) were eliminated from further analyses. Altogether, the recordings were very clean in all participants, but one (Patient 10). The amount of data rejected for this patient was almost 20%, whereas it was <1% for all other patients. Power spectral densities (PSD) were then calculated using a Short-Time Fourier Transform, with a 0.25s Hamming window and 60% overlap, and finally log powered.
We propose a model that describes the LFP spectrum and allows the extraction of two components of the PSD: the broadband activity and the beta oscillatory activity. Similar approaches have been used to model the PSDs in EEG 2,3 . The different components of the model are illustrated in Figure   1c (main text), in which the LFP PSDs are modeled as: where g( ) is the modeled power spectrum at frequency .
The function 1 ( ) models the broadband activity function at frequency with a power law function; 1 , 2 , 3 . Parameter 1 determines the constant offset of the PSD. Parameter 3 represents the speed of decay of the 1/f function, and can be visualized as the slope at f=1Hz.
Finally, parameter 2 allows to change the rotation axis of the 1/ 3 function. If 2 =1, the rotation axis of 1/ 3 will always occur at 1 µV 2 at f=1Hz, no matter the value of 3 . The speed of decay 3 defines how the PSD rotates around this axis. Parameter 3 allows to shift the rotation axis to a different value than 1 µV 2 at f=1Hz, which may allow to model more realistically neurophysiological phenomena 3 .
The function 2 ( ) models the beta activity function at frequency with a Gaussian probability density function, with amplitude 4 , mean frequency 5 , and standard deviation 6 . Fig. 1d (main text) exemplifies the impact of changing each of the parameters on the modeled broadband (top row) or beta activity (bottom row), respectively. 1−6 parameters were estimated from the entire 6-min recording by an interior-point optimization algorithm, with the objective function defined as the root mean square error between the modeled and actual PSDs: where frequencies ranged between 8-90 Hz, except for frequencies between 48-52Hz which were not taken into account for the modeling due to the use of a notch filter to remove power line noise. We have ignored frequencies below 8Hz in order to avoid low frequency artifacts 1 and tremor-associated activity 4   has the same beta amplitude as the light green curve, but the level of broadband activity is larger, resulting in an overall larger total amplitude. Finally, the dark green curve has the same overall amplitude as the neutral green curve, but the beta amplitude is actually much larger. These examples emphasize the need to model these physiological processes independently.
Although, we cannot rule out that the broadband activity might actually be a combination of device noise and true broadband activity, the same equipment was used for all the recordings performed.
As such, potential impedance differences should average out at the group level (patients and hemispheres), and results obtain should underlie physiological processes.
Emerging evidence has also shown the existence of two distinct beta sub-bands carrying pathological information 5 : the low-beta (13-20 Hz) and the high-beta (20-35 Hz) bands. As such, we also evaluated the feasibility to model the PSDs with two beta bumps. To this end, the PSD was modeled as the sum of two Gaussian function and a power law function, and therefore estimated nine parameters instead of six parameters in the case of a single beta bump.

Statistical analysis
To validate our model, we first evaluated the goodness-of-fit between the modeled and actual PSDs using the coefficient of determination ( 2 ). We then quantified the relationship between the model parameters 1−6 and the clinical scores using Spearman's correlation coefficient (rho), pooling patients and hemispheres (N=26). The significance level was computed using exact permutation tests and corrected for multiple comparisons using false discovery rate (FDR) 6 .
We then compared differences in correlation with clinical scores, when the beta oscillatory activity was estimated with and without broadband activity. The estimated beta activity without broadband To validate our approach, we performed the same analysis using the actual PSDs, and compared the results with those obtained using the modeled PSDs. For this, the beta peak in the actual PSDs for each patient and hemisphere was determined by an external expert blinded to the study. Then, we evaluated the correlation of the beta amplitude in the actual PSDs with motor symptoms, with and without the modeled broadband activity.
Finally, we investigated the relationship between the broadband activity and motor symptoms and age, as previous works have shown that age correlates with broadband activity within specific frequencies 7 . As such, we evaluated for each patient and hemisphere (N=26) the broadband activity function (function 1 , see Figure 1c in the main text) within the frequencies ∈ [0,200] Hz, based on the parameters estimated from the model. Then, for each evaluated frequency , we correlated comparisons using the false discovery rate method (FDR) 6 . For the correlation of the broadband activity function with age, we associated the same age to both hemispheres for a given patient.

Code availability
Data analyses were conducted in Matlab using scripts available from the corresponding author upon reasonable request.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Modeling of the power spectrum
Goodness-of-fit ( 2 ) between the modeled and actual PSDs ranged between 0.95-0.99, and thus our mathematical model accurately fitted the power spectrum (see Supplementary Figure 1). While Parkinson's disease has often been associated with two beta frequencies (low beta (8-20Hz) and high beta (21-35Hz)) 5 , our rationale for performing the analysis on a single bump was two-fold: first, there were no significant differences across conditions in the adjusted R 2 -which takes into account the number of parameters as a regularization method (mean R 2 oneBeta=0.98 and mean R 2 twoBeta=0.98, p = 0.16; signed rank test); and second, we found no significant differences across conditions in the Akaike information criteria, and it was lower on average for the case of one beta bump (mean AIConeBeta= 496 and mean AICtwoBeta = 525; p = 0.43; signed rank test). These results suggest that power spectral densities can be modeled using customizable parameters depending on the needs. Despite this, and to avoid an overfitting effect when modeling two bumps, we limited our investigation to modeling beta activity using only a single beta bump.

Correlation between model parameters and motor symptoms
In order to identify how the LFP power relates to motor symptoms, we investigated the relationship between the different parameters of the modeled power spectrum and the UPDRS III clinical scores. The fitted amplitude of the beta activity (k4) correlated with both bradykinesia (rho=0.55, p=0.006) and rigidity (rho=0.69, p=0.000), but not with tremor (rho=-0.07, p=0.743; Figure 2 main text). In addition, the beta amplitude correlated with bradykinesia plus rigidity (rho=0.68, p=0.000) and all three clinical scores summed together (rho=0.50, p=0.011). No other parameters presented significant correlations (p>0.1; Supplementary Table 3). These findings agree with previous electrophysiological studies showing that pathological beta activity was associated with bradykinesia and rigidity, but not tremor 8 .
Similar to the results obtained with our model, the correlations between motor symptoms and beta amplitude using actual PSDs were significant only when removing the broadband activity (FDR corrected p<0.05 for Bradykinesia, Rigidity, B+R and B+R+T); but they were not significant when with and without broadband activity was significant for Rigidity (Hotelling's t-test, FDR corrected p=0.02) and B+R (FDR corrected p=0.03).
To further evaluate the validity of our model, we compared the beta amplitude and beta frequency estimated by the model with those obtained using the actual PSDs. Results showed a very strong correlation between the beta power in the real data and the modeled beta peak (Supplementary . Similarly, the correlation between the modeled beta frequency and the real beta frequency was also very high (r=0.62, p=0.0006), further highlighting the accuracy of the model. Figure 1. Model fits to individual datasets. Model fits (orange = broadband activity, blue = beta activity and green = sum of broadband and beta activity) for each patient and hemisphere, together with the actual power spectral densities (black lines). For each subject and hemisphere, the bradykinesia, rigidity and tremor scores are also reported, as well as the individual R 2 goodness of fit. Note that the fitting was not performed (nor shown) within the 48-52Hz due to the use of a notch filter.