Cortex-wide topography of 1/f-exponent in Parkinson’s disease

Parkinson’s disease (PD) is a progressive and debilitating brain disorder. Besides the characteristic movement-related symptoms, the disease also causes decline in sensory and cognitive processing. The extent of symptoms and brain-wide projections of neuromodulators such as dopamine suggest that many brain regions are simultaneously affected in PD. To characterise brain-wide disease-related changes in neuronal function, we analysed resting state magnetoencephalogram (MEG) from two groups: PD patients and healthy controls. Besides standard spectral analysis, we quantified the aperiodic components (κ, λ) of the neural activity by fitting a power law κ/fλ – f is the frequency, κ and λ are the fitting parameters—to the MEG power spectrum and studied its relationship with age and Unified Parkinson’s Disease Rating Scale (UPDRS). Consistent with previous results, the most significant spectral changes were observed in the high theta/low-alpha band (7–10 Hz) in all brain regions. Furthermore, analysis of the aperiodic part of the spectrum showed that in all but frontal regions λ was significantly larger in PD patients than in control subjects. Our results indicate that PD is associated with significant changes in aperiodic activity across the whole neocortex. Surprisingly, even early sensory areas showed a significantly larger λ in patients than in healthy controls. Moreover, λ was not affected by the Levodopa medication. Finally, λ was positively correlated with patient age but not with UPDRS-III. Because λ is closely associated with excitation-inhibition balance, our results propose new hypotheses about neural correlates of PD in cortical networks.


INTRODUCTION
In Parkinson's disease (PD), the progressive loss of the dopaminergic cells not only depletes the neuromodulator dopamine but also alters the dynamics of other key neuromodulators such as serotonin, noradrenaline and acetylcholine (see the review by McGregor and Nelson 1 ). Given the widespread prevalence of neuromodulators in the neocortex, it is expected that the neocortical neural circuits dynamics and function would also be affected in PD. That is, the signature of PD-related dysfunction should also be visible in the neuronal activity recorded from neocortical regions.
Consistent with this, analyses of electroencephalogram (EEG) and magnetoencephalogram (MEG) have revealed several changes in the population activity of different neocortical regions (see reviews by Geraedts et al. 2 and Boon et al. 3 ). While the results are diverse given the heterogeneity of the patients, it is commonly observed that the low frequency-delta to low-alpha bandpower increases whereas high-alpha to gamma band power decreases 4 .
This kind of spectral slowing was shown to be correlated with motor and cognitive symptoms 5 . Such spectral slowing has been observed in the earliest stages of the disease (e.g., in the posterior cortex regions 6 ), hence demonstrating that it is not an effect of dopamine medication. Moreover, dopamine replacement therapy hardly reverses the spectral slowing, especially in the more advanced PD patients 6 . Compared to this (power) spectral slowing, the frequency peaks were much less studied. Results on the latter usually hold for global quantities such as the average of cortical oscillations (e.g., see Soikkeli et al. 7 ), or the whole spectrum dominant peak (e.g., see Geraedts et al. 2 for a review), which both decrease in PD patients. Here we study the change in frequency peaks which is somewhat similar as the (power) spectral slowing as frequency peaks follow high power parts of the power spectral density (PSD). This analysis is a way to test if spectral slowing is also seen in our data.
Spectral power alterations are also associated with changes in functional connectivity (FC) in PD patients. Throughout the disease, the low-alpha band FC decreases after an initial increase 8,9 . Increase in beta band synchrony is also commonly observed in the basal ganglia (BG) as well as in the cortico-basal ganglia loops (see the review by Hammond et al. 10 ).
Thus, previous work on analysis of EEG/MEG has been largely focused on oscillatory activity in different frequency bands. Besides oscillations, the aperiodic part of the population activity can also be informative about the underlying network dysfunction and relative excitation-inhibition balance 11 . To the best of our knowledge, the aperiodic part of the EEG/MEG activity in PD was characterised in three studies. First, Vinding et al. 12 reported steeper power-law decay of MEG power in sensorimotor regions of PD patients. In a more recent study, Wiesman et al. 13 estimated the aperiodic activity on four different frequency bands showing a neurophysiological slowing-decrease through frequency of the λ deviation from healthy control (HC). Finally, Wang et al. 14 studied low-spatial resolution EEG from PD-ON and OFF medication showing that there is a significant increase of the offset and exponent power-law parameters from OFF to ON medication in some of the sensors. However, thus far it has remained unclear how the spatial distribution of aperiodic component of the population activity is altered by chronic dopamine depletion and dopamine replacement therapy in human patients. 1 Therefore, we studied the spatial distribution of spectral peaks and aperiodic component of the activity from neural populations measured with MEG in PD patients in their ON and OFF medication states. To this end, MEG acquired using 306 sensors were pre-processed to obtain 44 sources' activity distributed all over the neocortex according to the HCP-MMP1 atlas 15 . We found that indeed there is a slowing in the MEG spectrum even when the spectral peaks were identified without classical frequency band definition. Analysis of the power-law exponent (λ) of MEG spectrum revealed a significant increase in the λ in sensory and motor regions in PD patients compared to the healthy controls. In fact, λ showed a spatial positive gradient from anterior to posterior brain regions in PD patients. Surprisingly, the frontal regions which receive most of the dopaminergic projections did not show a significant difference in λ. Levodopa administration did not affect the spatial distribution of λ. However, λ changes were correlated to the patients' age but not to their Movement Disorder Society's UPDRS-III score (clinical rating of motor symptoms 16 ). That is, neocortical activity dynamics is more vulnerable to chronic dopamine changes in older PD patients than in younger patients. Moreover, our results reveal aspects of neocortical population activity which are not affected by Levodopa therapy. Because 1/fexponent (λ) can be linked to excitation-inhibition (EI) balance, our analysis suggests new testable hypotheses about the PD-related changes in the neocortex.

RESULTS
To characterise how cortical activity is changed in human PD patients, we analysed the MEG signals acquired during resting state. To this end, we quantified both the oscillatory peaks and the slope of the power spectrum for MEG sources corresponding to 44 different brain regions, and then we estimated the dependence of the latter on age and UPDRS-III, see "Methods" for details.
Slowing of spectral peaks in PD Analyses of MEG and EEG from PD patients suggest that frequency associated with peak power in the theta and alpha bands is reduced in human PD patients when compared to healthy controls 3 . To test whether this is also the case in our data, we searched for Gaussian peaks in the spectrum of MEG signals without explicitly defining the frequency bands (see "Methods").
By estimating the distribution of all the frequency peaks observed across all the brain regions, we found that indeed in PD there is a general slowing of spectral peaks (PD-OFF vs HC, P value < 10 −20 with Kolmogorov-Smirnov test, see Fig. 1a). Notably, Levodopa medication did not improve this slowing, in fact, if at all it seemed to worsen the frequency peaks slowing (PD-OFF vs PD-ON, P value < 0.017 with Kolmogorov-Smirnov test).
Next, we sorted the frequency peaks into classical frequency bands (theta = 4-8 Hz, alpha = 8-13 Hz, beta = 13-30 Hz, and gamma = 30-45 Hz). Thus, we obtained the proportion of peak frequency (p band;b G ) for each band, MEG source b and group G. Given the frequency peaks slowing there were not enough peaks in the theta and gamma bands to compare the HC and PD groups (Fig. 1a, c, d). This explains the missing dots for some brain regions in Fig. 1c, d. Therefore, we restricted our statistical analysis to alpha and beta bands as it is not possible to perform statistical analysis with the few points we have in theta and gamma bands.
Even though there was a reduction in alpha band peak power, the difference did not reach a statistical significance level for none of the brain regions (Fig. 1b, left). In the beta band also we found a widespread decrease in the mean peak frequency among most of the patients (Fig. 1b, right). This suggests that frequency peaks slowing does not affect the alpha and beta bands locally. However, when we pool data across all the brain regions, frequency peaks slowing can be observed in all bands. The main local effect (i.e., brain region-specific) of frequency peaks slowing is the reduction in the power of gamma-band oscillations and increase in the power of theta bands oscillations. In particular, gamma-band peaks were almost completely missing from the frontal regions in PD patients (Fig. 1d). By contrast, HC seem to be lacking theta band peaks in the central-temporal and posterior brain regions (Fig. 1c).
The lack of significance might be due to the low number of data points which makes it very hard for changes in frequency peaks distribution to be significant. But it is also possible that each patient's brain compensates in its own way, and therefore we do not find statistically significant changes at the group level. Because the periodic component of the data was not informative about the disease, we performed the analysis on the aperiodic component in order to get more insights on how these changes are spatially and temporally distributed over the cortex.
Neocortex-wide change in the 1/f-exponent distribution in PD Next, we focused on the aperiodic component of the neural activity and analysed how the spectral power decreased as a function of frequency. Using the FOOOF algorithm, we fitted a power-law function to the PSD (see "Methods") and estimated the exponent λ for each MEG source and subject. Across all the data, λ spanned a relatively wide range 0.1-2.0. However, temporal variation of λ for individual brain regions in both healthy controls and PD patients was smaller than across subject variance ( Fig. 2a, b). We found that in healthy controls, λ was smaller in sensory regions than in the cognitive regions (Fig. 2c, left). However, this apparent gradient of λ from frontal to posterior regions is not statistically significant.
By contrast, in PD patients (both OFF and ON medication) λ was larger in sensory regions than in the cognitive regions (Fig. 2c, middle, right). Moreover, λ showed a clear frontal to posterior gradient (Spearman correlation (r s ), P value < 10 −6 and r s ðλ b G ; y b Þ < À0:65 where y b is the y-coordinate of the brain region b, G ∈ {PD OFF, PD ON} and λ stands for the time average of λ).
From Fig. 2c, it is clear that there are differences in the spatial distribution of λ in PD patients and HC. To quantify this difference, we performed a brain region-wise comparison of λ in HC and PD patients OFF medication (Fig. 2d, right). We found that PD patients (OFF medication) have larger λ in auditory, visual, somato-senssory and motor regions than in HC (for significant figures, please see Fig. 2d, right). The same landscape of differences in λ was observed when we compared PD patients ON medication with HC ( Fig. 2d, middle). Interestingly, a comparison of PD patients in ON and OFF medication states did not reveal any significant differences across the whole neocortex (Fig. 2d, left). We also note that still in most brain regions (except frontal areas), the group variance among HC was much lower than among PD patients (edge colours of the dots in Fig. 2c). Finally, the most significant changes were observed in the left hemisphere, which is expected as in our cohort most patient's disease started in their right hand.
λ is not constant across time (Fig. 2a, b). To quantify these fluctuations over time, we measured the coefficient of variation (CV, see "Methods") of λ for each brain region in each subject. Small value of CV of λ indicates temporal stability of λ. In general λ was stable over time (small CV). However, λ was more variable in frontal regions in both PD patients and HC ( Supplementary  Fig. 4c). A comparison of CV of PD patients and HC revealed that for most brain regions λ was less variable in PD patients than in HC, particularly in the sensory regions where the CV difference was statistically significant ( Supplementary Fig. 4d, right).
Moreover, even if there is no significant difference in this second-order temporal statistics between PD patients ON and OFF medication, Levodopa seems to decrease temporal fluctuations in Fig. 1 Spectral slowing, a peak frequency perspective. a Distribution of peak frequencies in all brain regions of all individuals from the different groups. b Distribution of frequency peaks in alpha and beta bands. The colour inside the dots is the difference (PD-OFF -HC ses1) of the mean peak frequency of each brain region in alpha (left) and beta (right). The edge colour of the dots is the normalised difference (PD-OFF -HC ses1) of proportion of individuals having at least one peak frequency in the given band. c, d Distribution of frequency peaks in theta (c) and gamma (d) bands. Dot colour indicates mean peak frequency within the band. Edge colour of the dots indicate the proportion of individuals having such peaks in the given group. some brain regions, especially frontal ones ( Supplementary  Fig. 4d, left).
Finally, as λ is relatively stable over time, one can estimate λ for the whole time series instead of small overlapping epochs. However, a single λ value for the whole time series obscures the temporal fluctuations in the aperiodic component which may be informative about the disease-specific brain dynamics. Indeed, we observed that fluctuations are weaker in PD patients even though they have higher aperiodic components ( Supplementary Fig. 4c, d). This is counterintuitive as we might expect larger variations for higher λs.
The 1/f-exponent is positively correlated with age but not with UPDRS-III in PD In the above, we have shown how λ was altered in PD patients across the whole neocortex. The question arises: does this difference in λ also correlate with the PD severity? λ could be affected by subject age and UPDRS-III score. Therefore, we estimated partial correlations between λ and UPDRS-III given age, and age given UPDRS-III for each brain region separately (see "Methods").
We found that a negative partial correlation between λ and UPDRS-III scores but these correlations do not reach a statistical  Table  1 for exact P values and Supplementary Fig. 6 for Spearman and Pearson correlations). By contrast, we did not find any significant correlation between λ and age in healthy controls. Despite a decline in the dopaminergic neuron function due to normal ageing 17 , these "normal" changes do not seem sufficient to affect the aperiodic component of cortical population activity, unlike in PD.
It is somewhat counterintuitive that despite significant PDrelated changes in λ, it was not correlated with UPDRS-III, even in the motor regions of the cortex. This result may be explained by a lack of sensitivity of the data or by the fact that λ is not by itself related to motor symptoms. The fact that the increase in λ was positively correlated with age suggests that persistent dopamine depletion may be more detrimental for brain networks that are already undergoing age-related deterioration like the sensory areas. For instance, vision is impaired in PD, especially in patients with cognitive decline 18,19 . Alternatively, brain networks in younger patients may be more malleable to compensate for the loss of dopamine than in older patients.

DISCUSSION
Given chronic change in the neuromodulator dopamine, we expect that excitation and inhibition (EI) balance may be altered in the neocortical regions. To test this we analysed the aperiodic activity of MEG data which has been inversely related to EI balance in several previous studies [20][21][22][23][24][25][26][27][28][29] . We found that the 1/f-exponent (λ) was higher in PD patients than in HC in most brain regions but not in frontal regions. The most significant changes occurred within sensory and motor regions. Moreover, we found that λ showed a spatial gradient from posterior to anterior brain regions in PD patients. Although not significant, in HC the gradient of λ was reversed compared to that in the PD patients. When examining the fluctuations of λ over time, we found that λ was more variable in the frontal regions. Moreover, λ fluctuated more in HC controls than in PD patients. This was not expected as λ was globally larger in PD. Surprisingly, Levodopa medication had a little effect on the topography of λ. Finally, we computed different correlation measures between λ and both the age and UPDRS-III score. We did not find any correlations in the HC group either with age or UPDRS-III. Although not significant, global negative correlations with UPDRS-III were observed in PD patients. On the other hand, λ was positively correlated with age in PD patients in both ON and OFF medication states.
Previously, Vinding et al. 12 analysed the same data for 1/f slope (λ) but that analysis was restricted to the sensorimotor association regions. Consistent with our results, they showed that the exponent is steeper in the sensorimotor cortex of PD patients. Still aligned with our results, they also found that λ was positively correlated with age in PD patients but not in HC and the correlation between λ and UPDRS-III was not significant. We build on that study and now show the topography of λ over the whole neocortex in both healthy controls and PD patients. Our analysis shows that λ changes in early sensory regions are also correlated with age but only in PD patients.
Recently, Wang et al. 14 analysed the 1/f slope of EEG activity. Our results and methods differ from theirs in several ways. First, Wang et al. 14 acquired brain activity using 32 EEG sensors, whereas we have a much denser spatial sampling of brain activity as we used 306 MEG sensors. Moreover, we performed a source reconstruction whereas their analysis was at the sensor level. Concerning the topography of λ's topography, we found that λ decreased (increased) from sensory to cognitive regions in PD patients (healthy controls). By contrast, Wang et al. 14 have reported similar topography in both healthy and PD patients. In their study, λ was high in the central regions with a decrease in all directions from there on. This difference might be due to the fact that we used a brain atlas to project the MEG sensors to sources as well as the difference in the number of EEG/MEG channels used in two studies. Nevertheless, this is a conspicuous difference that should be investigated in a study where both MEG and EEG are acquired from the same patients. In contrast to Wang et al. 14 , we did not find significant differences between PD-ON and PD-OFF states. This could be due to the difference in recording protocols in the two studies: in our case, dopamine effects were acute in the sense that patients were OFF medication for at least 12 hours but then took the medication only one hour before the ON medication state was recorded. Unlike ref. 14 , we have analysed partial correlation between UPDRS-III/age and each MEG source. This separation of brain regions and partial correlation revealed that λ was correlated to age but only in PD patients.
There are at least three possible interpretations of the changes in λ. Population signals such as MEG are generated by dipoles created by transmembrane currents 30 . The transmembrane currents are generated due to synaptic inputs impinging on the neurons. Therefore, the frequency spectrum of the MEG (and also EEG, LFP) reflects the time constants of the synaptic inputs. Typically, inhibitory (GABAergic) synaptic transients have a longer time constant than excitatory (glutamatergic) current. Therefore, increased λ may be an indicator of either increased GABAergic or reduced glutamatergic currents. That is, λ is a proxy of relative excitation-inhibition (EI) balance 11 .
However, even glutamatergic synaptic inputs due to NMDA receptors can also be very slow. Therefore, it is also possible that increased λ may indicate a relative increase in the NMDA type synaptic currents. In the context of our results, we found that λ in PD is increased in the sensory areas. These areas usually have relatively small fractions of NMDA receptors 31 .
Finally, it is also possible that the change in λ is simply a reflection of changes in the dynamics of local and external inputs to the network. That is, if the network is driven by slowly fluctuating inputs it will also reflect in slower fluctuations and Fig. 4 Cortex-wide projection of the centre of mass of the 44 brain regions investigated here using MEG. The brain regions were extracted using the HCP MM1 atlas. therefore larger λ. If this is the case then we would also expect a change in the time scales of spiking activity.
Each of these possibilities suggest a different but testable hypothesis. If the slope of the MEG reflects change in the relative fraction of excitatory and inhibitory currents, then it is tempting to hypothesise that in PD there is an increase in inhibition in the sensory region. Increase in relative fraction of NMDA will also account for an altered excitation-inhibition balance on longer time scales. Techniques such as MRS could be used for a non-invasive estimate of the relative fraction of AMPA, NMDA and GABA in order to test this hypothesis. Ideally, these hypotheses should be tested in animal models with in vivo measurements of excitatory and inhibitory currents.
If changes in the λ reflect a change in the time scales of fluctuations in the input then, we will need to explain how in a presynaptic network the spectrum of neuronal activity could change. For that we again revert to the hypothesis of a change in excitation-inhibition balance. However, before following this line of reasoning, we need to estimate whether there are significant changes in the spiking activity of a given brain region where λ has changed. This suggests that spiking activity from early sensory regions should also be recorded in animal models of PD.
As we discussed earlier, it is tempting to relate the slope of the frequency spectrum to EI balance. Origin of pathological activity in the basal ganglia during PD, especially the beta band oscillations is closely related to changes in the EI balance in STN and GPe regions of the basal ganglia 32,33 . Local field potential recorded from the STN or GPe during DBS surgery can be used to estimate relative EI balance at different locations in the STN in terms of λ s. Such an estimate of relative EI balance could guide stimulation electrode placement.
A lack of a clear correlation between λs and UPDRS-III suggests that it may not be useful as a clinical bio-marker. The largest changes in λ were observed in the sensory areas. Therefore, it is reasonable to hypothesise that λ should be related to the severity dysfunction of sensory processing or sensorimotor integration. But it should be emphasized that UPDRS-III does not measure sensory dysfunction. This could be one of the reasons why λ is not related to the severity of PD motor symptoms in our study.
However, the correlation between λs and age in the PD group suggests that λ is indeed altered in PD. In particular, we observe a positive correlation in PD in most brain regions except the frontal regions. This result may suggest that dopamine depletion has an important impact on the sensory and motor "ageing". In particular, this effect seems to be non-linear and non-monotonic with respect to age. Indeed, we used different dependence measures to tease apart the different possible dependence types: linear (Pearson), monotonic (Spearman), no assumption (distance correlation). Our results show that when there is a dependence, it seems to be mainly non-linear and non-monotonic (significant and larger distance correlation than other dependence measures). This is somehow expected as distance correlation theoretically includes any kind of dependence but here it is much more significant in the age dependence case.
A positive correlation between λ and age in PD patients is a curious observation. Previous work from several groups have shown a negative correlation between λ and age in the neocortex [34][35][36] . However, in over 60 years old HC, λ from somatosensory regions may positively correlate 37 . In contrast to these findings, in our data we did not observe any significant correlation between λ and age in HC. As age is highly (positively) correlated to disease duration, PD thus seems to affect all (except frontal parts) necortex's λs by increasing them over age/disease duration.
Levodopa effects are fast as shown by the UPDRS-III score improvements only one hour after the drug administration. We found that in our data Levodopa influenced the frequency peaks slowing. However, the effect of Levodopa on λ was not observed in the neocortex. To the best of our knowledge, there is no consensus on the cortical effects of Levodopa 12 . In our study, a possible reason could be the short time (1 hour) between taking the medication and recording of MEG. A recent study in ref. 14 reported changes in the λ estimated from EEG. In that study, PD patient ON medication took their medication dose as usual, in the morning before the measurement. However, Wang et al. 14 also reported that Levodopa did not improve the aperiodic part (i.e., λ increased) of the neocortex activity spectrum. In addition, the total washout of Levodopa may take days 38 . The latter could explain the similar results we obtained for ON and OFF medication.
Here we used a commonly used range of frequencies (1-45 Hz) for our analysis. This range is relatively small. Indeed the broader the frequency range, the better the fit should be. However, for MEG it is not as easy to take the largest band possible because muscle artifacts become more prominent in high frequencies. Therefore, local field potential or ECoG may be more suited for such an analysis over broad frequency ranges.
Next, we have used a specific brain atlas and it is not clear how the topography of λ may change when using a different brain atlas. In this regard, adding information from the sensor space may help matching the different atlases' results.
Finally, to better understand the functional implications of λ changes, a more detailed correlation analysis is needed that takes into account the UPDRS-III sub-scale 16 . Furthermore, future studies should also take into account the non-motor symptoms by adding cognitive, depression and anxiety scores in our analysis. However, as the latter are quite difficult to evaluate, we might need a larger cohort than ours to interpret correctly the results obtained from their analysis.
Our findings also raise the question why frontal regions are more protected than sensory and motor regions even though dopaminergic projections in the neocortex are primarily restricted to the frontal regions 39 . The lack of prefrontal changes between PD and HC may be due to dopamine pathways. Indeed, the main brain area producing dopamine and projecting to the cortex (in particular prefrontal cortex) is the Ventral Tegmental Area, which is altered after the substantia nigra compacta in PD 40 . Hence, changes in the frontal regions may take longer to manifest.
In addition, because there is such a big change in the λ value in sensory regions one would expect deficits in the sensory representations. It is well established that PD patients have olfactory 41 , proprioceptive 42 and cross-modal sensory fusion deficits 43 . However, our work suggest changes in other sensory modalities such as vision and audition.
It is common to measure functional connectivity between brain regions in a frequency-dependent manner. Usually these estimates are based on filtered time series. Here we found that λ varies over time. Therefore, we can ask whether λ variations across brain regions are correlated or not in PD and HCs. Recent work suggests that functional connectivity based on the component of aperiodic activity may be more robust 44 .
Overall, we show that the aperiodic activity, which usually has been considered as noise, gives new insights in PD and deserves more attention when analysing any neural field potentials like ECoG, LFP, EEG and MEG. Finally, previous studies showed frequency peaks slowing in relation to cognitive decline rather than the motor symptoms 5 . It could then be that the change in λ is a more general expression of neurodegeneration than only the dopamine-affected systems.

MEG data and its pre-processing
Here we analysed the resting state (eyes opened) MEG recorded from 17 PD patients (age 41-85; five female) and 20 age-matched healthy controls (HC; age 54-76; eight female). The data were acquired at the Swedish National Facility for MEG (NATMEG, https://natmeg.se/) using the Elekta Neuromag TRIUX 306-channel MEG system. The study was approved by the regional ethics committee (Etikpöingsnöamden Stockholm, DNR: 2016/911-31/1) and followed the Declaration of Helsinki. All participants gave written informed consent before participating.
For more details of the experimental protocol, please refer to ref. 45 study. Briefly, the MEG was acquired at a sampling rate of 1 kHz with an online 0.1 Hz high-pass filter and 330 Hz low-pass filter. Each subject was recorded twice. MEG from PD patients was recorded in OFF and ON medication states. In the OFF medication condition, PD patients were off their dopamine replacement medications (Levodopa) for at least 12 h. After the first recording session, patients took their medication and MEG was recorded for the second time one hour after the medication intake. HC subjects were also recorded twice at an interval of 1 h. Neurological and motor function MDS-UPDRS-III 16 (shortened to UPDRS-III or only UPDRS in the following) tests were performed on all participants between the two sessions. PD patients had a second motor test after the second session (ON medication). The patients included in this study represent a relatively homogenous non-demented PD sample in early to moderate disease stages who showed a very good response to dopaminergic medication (supported by a mean reduction of UPDRS-III by~50%). A summary of these results is given in Supplementary Table 2 as well as disease duration and Levodopa dose.
Each recording epoch was of 8-min duration. Data were downsampled to a sampling frequency of 200 Hz. Downsampling was done to speed up the analysis process without altering the results. Further, we low-pass filtered the MEG signals (cutoff frequency 45 Hz). This was done because MEG signals do not have a high signal-to-noise ratio at high frequencies. We pre-processed the data to remove non-neural noise and eye movement artifacts using ICA and data from an electrooculogram or electrocardiogram. We then used the dynamic statistical parametric mapping (dSPM 46 ) for source reconstruction (implemented in MNE-Python 47 ), followed by labelling using HCP MM1 atlas 15 which resulted in 44 brain regions (BR) signals as shown in Fig. 4.

Spectral analysis
The spectral analysis was done with the FOOOF toolbox 48 . We estimated two main parameters from the spectrum: the frequency of spectral peaks (f peak ) and the aperiodic components (see below). f peak was measured for the whole recording, whereas the aperiodic components were extracted for short epoch of signal to obtain its temporal dynamics.
Extraction of the frequency peaks. For each brain region of each patient we estimate the power spectral density (PSD)using Welch's method 49 implemented in SciPy 50 with default parameters (average = "mean" and window = "hann")averaging over 5 s segments with 50% overlap over the whole signal. The periodic part of the spectrum was extracted using the FOOOF method proposed in ref. 48 .
FOOOF optimally fits the PSD with a function composed of the sum of an aperiodic part and a periodic part. The aperiodic part is modelled as κ f λ where f > 0 Hz is the frequency, κ is the offset of the PSD, and λ defines how the PSD decays as a function of f. The periodic part is modelled as a sum of P weighted Gaussian distributions whose mean, standard deviation and weight, respectively, indicate the peak frequency, width and height of the spectral power bump of the peak. We manually specified the frequency band (1-45 Hz) on which the fitting was done. To identify oscillatory peaks in the PSD using FOOOF, we need to provide the maximum number of peaks P, their minimum height and bandwidth. We searched for P = 4 peaks with a minimum peak height of 10 0.2 V 2 /Hz and peak width within [1,10] Hz (see Fig. 5c). Thus, we obtained f 1 peak ; ; f 4 peak ; for each channel and subject. For Fig. 1a, we pooled all f 1 peak ; ; f 4 peak ; for all channels and subjects in each group.
From the proportion of individuals having such peaks within each group, we defined a normalised difference between groups as ðbÞ is the proportion of individual in the group G ∈ {PD OFF, HC} having at least one peak frequency in the frequency band band and the brain region b (see Fig. 1b). These proportions, p band G ðbÞ, for the theta and gamma frequency bands are illustrated in Fig. 1c and d, respectively.
Calibration of λ estimate. A priori it is not clear how much error is expected in the estimation of λ given a specific λ and epoch size.
We have fixed the epoch size to 50 sec, which is large enough to give us a good estimate of the PSD. Still we need to determine the appropriate segment size (W) for Welch's algorithm. To address this issue, we resorted to a numerical simulation approach. First, we created a PSD of the form κ f λ true , and assigned random phase (drawn from a uniform distribution U ½Àπ;π ) to each frequency. Next, we used the inverse Fourier transform to reconstruct a signal of desired length ( Supplementary Fig. 1a). Finally, we extracted λ sim for different values of W ( Supplementary Fig. 1d-f). We systematically varied λ true and W. We found that the error in the estimate of λ sim depends on both the λ true and W ( Supplementary  Fig. 1e, f). Based on these numerical experiments, we chose W = 2 sec. which gave us the best compromise between the errors in mean and standard deviation of the λ in a range that we have observed in our data.
Temporal dynamics of the slope of the MEG power spectrum. To characterise the aperiodic component of the MEG activity, we focused on the way the power of the MEG signals decayed as a function of the frequency. To this end, we first estimated the time-resolved spectrum (spectrogram) of each MEG source (epoch size = 50 s, overlap 40 s, see Fig. 5a). The PSD for each epoch was evaluated using Welch's method averaging over 2-s segments with 50% overlap (default parameters: average = "mean" and window = "hann"). Next, we fitted a power law to each PSD using the FOOOF algorithm. In the aperiodic activity analysis, we used the frequency band [1,45] Hz, a minimum peak height of 10 0.2 V 2 /Hz, a peak width within [2,10] Hz and a maximum of 4 peaks (see Fig. 5c for an example of such a FOOOF fitting). The peak width lower bound is larger than the one used in the frequency peaks analysis because the frequency resolution is lower here. It is important to have a lower bound significantly larger than the frequency resolution to avoid confounding frequency peaks due to noise. Here we used 4 (resp. 5) times the frequency resolution for the aperiodic (resp. frequency peaks) analysis. Finally, we obtained the temporal and spatial dynamics of λ, κ and P by applying this method to each epoch of each source (see Fig. 5d red trace). We checked the accuracy of the fit by computing its r 2 error (for all the data combined; mean: 0.94, standard deviation: 0.04).

Variability of λ
For each brain region, we assume that λ depends only on two parameters, time and group (PD patient ON or OFF medication and HC). So, for a brain region b and a subject i; λ b i ðtÞ is the typical dynamics of the 1/f-exponent. To estimate the variability of λ b i ðtÞ as a function of time in a given individual and BR, we calculated the coefficient of variation (CV) as: where σ λ b i and λ b i are, respectively, the standard deviation and the mean of λ b i ðtÞ.
Relationship between λ and subject age and UPDRS-III scores λ could be related to the subject age and disease severity (UPDRS-III score) in both linear or non-linear manner. Therefore, we estimated both linear and non-linear measures of dependence between λ and x ∈ {age, UPDRS} using the following three descriptors: • Distance correlation measures the non-linear dependence between two variables: where using the notation || ⋅ || for the Euclidean norm, a jk ¼ jjλ j À λ k jj; b jk ¼ jjx j À x k jj; and 8c 2 fa; bg; c j ¼ 1 n Here λ i refers to the slope of the frequency spectrum of ith subject for a given brain region. Other variables and this definition are illustrated in Supplementary Fig. 2. For further details, please see the work by Székely et al. 51 .
• Pearson correlation measures the linear link between two variables: where cov is covariance and σ x is the standard deviation of the variable x.
• Spearman correlation measures the monotonic relationship between two variables: where R( ⋅ ) is the rank function.
To disentangle the dependence of λ on the age and UPDRS-III, we computed the partial correlations. To do so, we used partial_corr function of the pingouin 52 Python package to calculate the Pearson and Spearman partial correlations. Concerning the partial distance correlation, it is similar to the classical partial correlation. It is based on projections but now in a more complex space called the Hilbert space of U-centred matrices. We refer to the work by Székely and Rizzo 53 for a rigorous definition. We used the corresponding function, partial_distance_correlation, from the Python package dcor 53 . We performed a permutation test (200 000 permutations on the x variable) to estimate the P values.

Statistical tests
With the exception of the correlations among UPDRS-III, age and λ, we used the Kolmogorov-Smirnov test from the ks_2samp function implemented in SciPy 50 to determine the statistical significance of our results. This test captures more the deviations near the distribution centre than at its tails. Moreover, its power is greater when used in the one-tailed case. Therefore, we used the latter to test whether the cumulative distribution functions (CDFs) of one variable is greater or less in one group compared to another. In the test, the statistic is D þ ¼ sup u2R ½FðuÞ À GðuÞ where F and G are the CDFs to be compared 54 . For example, when comparing the λ in Fig. 2b, we tested whether the CDF of λ PDOFF was less than λ HC in the frontal regions and the opposite in other brain regions. When comparing the two HC groups-between sessions 1 and 2, see Supplementary Fig. 4bλ HCses1 and λ HCses2 are statistically closed (distributions across the group) so that we gathered both groups as one group (the HC group). To confirm this similarity, we estimated Pearson correlation between λ HCses1 and λ HCses2 , see Supplementary Fig. 3a. However, the CVs difference between the two sessions was of the same order than the difference between the CVs of HC and PD-OFF or PD-ON medication. Thus, we only compared CVs session-wise: PD-OFF and HC ses1, PD-ON and HC-ses2, see Supplementary Fig. 4c, d. Finally, in the frequency peaks analysis, because of the few number of frequency peaks in the theta and gamma bands, we could not do any meaningful statistical tests. Also, we only used session 1 for HCs as their frequency peaks are not reliable from one session to the other, see test-retest reliability of Supplementary Fig. 3b-f.
Concerning the statistical significance of correlation measures, we performed a permutation test (200,000 permutations on the x variable) to estimate the P values of the partial distance correlation. For Pearson and Spearman partial correlations, twosided P values were computed from the one-sample t test. The degrees of freedom of this test is N − 1 where N is the data size (N = 17 in PD groups and N = 20 in healthy groups).
Finally, we did not correct for multiple comparisons. We only used comparison-wise error rate for each specific brain area, because this study is an exploratory work on the possible changes of the 1/f-exponent dynamics in PD. The alpha level that was used to determine significance was a P value less than 0.01.

DATA AVAILABILITY
The MEG data cannot be made publicly available because of ethical permits. These data were collected under an ethical permit (DNR: 2016-19-31-1) and associated informed consents which states that "only participating researchers will access the data". For this reason, only researchers that are already included in the permit, or are added to the project by amendments to this permit, may access the data. In practice, this means that it is potentially possible, but not particularly easy, to add qualified researchers to the project, and thus allow access to the data. It should be noted that even if a new researcher is added to the project, all such access requires the data to stay where they are on the NatMEG servers at Karolinska Institute. For further questions, please contact daniel.lundqvist@ki.se.