Flexible brain dynamics underpins complex behaviours as observed in Parkinson’s disease

Rapid reconfigurations of brain activity support efficient neuronal communication and flexible behaviour. Suboptimal brain dynamics is associated to impaired adaptability, possibly leading to functional deficiencies. We hypothesize that impaired flexibility in brain activity can lead to motor and cognitive symptoms of Parkinson’s disease (PD). To test this hypothesis, we studied the ‘functional repertoire’—the number of distinct configurations of neural activity—using source-reconstructed magnetoencephalography in PD patients and controls. We found stereotyped brain dynamics and reduced flexibility in PD. The intensity of this reduction was proportional to symptoms severity, which can be explained by beta-band hyper-synchronization. Moreover, the basal ganglia were prominently involved in the abnormal patterns of brain activity. Our findings support the hypotheses that: symptoms in PD relate to impaired brain flexibility, this impairment preferentially involves the basal ganglia, and beta-band hypersynchronization is associated with reduced brain flexibility. These findings highlight the importance of extensive functional repertoires for correct behaviour.


Results
To test these hypotheses, we used source-reconstructed resting-state MEG data acquired from a cohort of 39 PD patients and 38 age-matched healthy controls (see "Materials and methods" for details). We analysed the size of the functional repertoire defined as the number of different unique avalanche configurations (see "Materials and methods") and estimated the contribution of individual brain regions to this measure of spontaneous cortical flexibility. Moreover, we correlated the individual flexibility of patients to motor and cognitive outcomes. Finally, we estimated the correlation between synchronization in the beta band (13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30) and the size of the functional repertoire. Data were acquired in two brief segments of eyes-closed resting-state conditions, each of 150 s duration. After artefact removal and segment rejection, the average duration of combined recordings was 131.05 s in the PD cohort and 131.46 s in controls. There was no significant difference in the average recording duration (Wilcoxon rank-sum test, p = 0.4539, difference of the average length = 0.41 s) or the number of ICA components removed in each group [KS-stats, p = 1.0000 and p = 0.9999 for electrocardiogram (ECG) and electrooculogram (EOG) components, respectively].
We first quantified the functional repertoire of these source-reconstructed resting-state MEG data, namely the total number of distinct avalanche configurations within each participant's data. A between-group contrast revealed that PD participants expressed a restricted functional repertoire, visiting a lower number of distinct patterns in comparable amounts of time (permutation test, p = 0.0088, see Fig. 2A,B). We seek to provide evidence of a putative mechanism in large-scale activity that could justify the restriction of the functional repertoire. Given the evidence that communication among brain areas can occur via synchronization, we reasoned that hyper synchronization across the network would indeed reduce its degrees of freedom, whereby the less flexible dynamics. Given that hypersynchronization has been consistently described in PD, we looked for a link between these two phenomena. Global synchronization was indeed higher in PD patients in the beta band-13-30 Hz (permutation testing, p < 0.001; Fig. 3A) and the size of the functional repertoire correlated negatively with the global synchronization in the beta band (r = − 0.4, p = 0.0070) (Fig. 3B). A negative correlation was found between the number of functional states visited and the clinical severity, as measured by the UPDRS-III, in the PD cohort (r(37) = − 0.31, p = 0.026) (Fig. 3C).
When relating the number of visited functional states with specific cognitive domains, evaluated by the MoCA subscores, linear positive relationships were present with the executive (r(37) = 0.34, p = 0.0352), attention (r(37) = 0.31, p = 0.049), language (r(37) = 0.32, p = 0.044), and naming domains (r(37) = 0.39, p = 0.015). We then compared the switch rate across all areas (i.e. the rate at which any region crosses the threshold), capturing the www.nature.com/scientificreports/ total number of activations regardless of their spatial configuration. Interestingly, PD patients showed a higher rate of switching (permutation test, p = 0.0047, Fig. 2C), implying that the decreased diversity of configurations in PD reflected more stereotyped activity, and not fewer activations per se. We moved on to study the heterogeneity of avalanches across members of the same cohort, and compared this heterogeneity between the groups. Permutation testing showed that the Hamming distance across avalanche patterns belonging to either healthy controls or PD differed significantly (permutation testing, p = 0.0048), where the avalanches in the PD patients were more homogenous in comparison to controls (Fig. 2D). This finding confirms that the ability to efficiently vary the activation patterns is impaired in PD patients. Finally, we estimated whether the distribution of the probability of each region to participate in the avalanche patterns differed between those patterns that were specific to a group vs patterns that were present in both groups. Using the Kolmogorov-Smirnov test, we observed that the two distributions were indeed statistically different (p = 0.012), implying a broad difference between regional involvement in avalanches in each group. Using bootstrapping to perform post-hoc analyses, we identified those regions that appeared more frequently in avalanche patterns that were unique to one group. The pallidum bilaterally (p = 0.03, p = 0.02, left and right, respectively), the left thalamus (p = 0.03), and the right putamen (p = 0.04) were more often involved in avalanches in the PD group than in the controls. Furthermore, the caudate bilaterally, the right thalamus and the left putamen showed differences that did not reach statistical significance (see Fig. 4).
We tested the robustness of our results to specific choices of the avalanche threshold and bin length by varying these variables across a moderate range of values and repeating the analyses. Specifically, we first used different binnings, ranging from 1 to 5 (i.e. each bin of size n results from n time points of the binarized time series). The results remained unchanged, with patients displaying a restricted functional repertoire for all the binnings explored (for the case of no binning (binning = 1) and binning equal to 5, p = 0.011 and p = 0.0068, respectively; see Supplementary Fig. 1). Furthermore, the avalanche threshold was modified, ranging from 2.5 to 3.5. For both cases the differences between the groups were confirmed (for z = 2.5 and 3.5, p = 0.0084 and p = 0.0082, respectively; see Supplementary Fig. 2). We also studied the impact of including cerebellar sources, and the difference in functional repertoire remained significant (p = 0.0042; see Supplementary Fig. 3). The analysis was also repeated for the classical frequency bands (i.e. delta, 0.4-4 Hz, theta 4-8 Hz, alpha 8-13 Hz, beta 13-30 Hz, An avalanche is a sequence of activations that starts when one or more regions are above threshold and ends when no region is above threshold anymore. The brains-plots for each moment (bin) of the avalanche show the areas above (yellow) and below (green) threshold (C). The areas that were above threshold at some time during the avalanche together form the avalanche pattern.  www.nature.com/scientificreports/ gamma [30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48]. In all the frequency bands, the difference in functional repertoire remained significant (p = 0.049, p = 0.020, p = 0.0065, p = 0.028, p = 0.026, from delta to gamma, respectively; Supplementary Fig. 3). These results suggest that the restriction of the functional repertoire in PD is a widespread phenomenon in terms of the range of frequencies involved. Finally, when repeating the analysis using the exact same amount of data for all subjects, the group difference in functional repertoire remained significant (p < 0.0001, for broadband; see Supplementary Fig. 4).

Discussion
Here we investigated time-resolved flexibility of brain activity in PD, utilizing avalanches and tools from statistical physics, adapted to MEG activity. The functional repertoire was formed by distinct patterns of activation during avalanches, and the size of the functional repertoire was used as an indicator of the ease of transitions between different brain functional states, and hence dynamic flexibility 7 . We found that PD patients exhibited a reduced functional repertoire compared to healthy controls. Moreover, the clinical disability was more pronounced in patients whose functional repertoire was more impoverished. Finally, we demonstrated the robustness of our results across a range of parameters and frequency bands, showing that the results did not depend on the specific details of data processing. Our results showed that PD patients display a restricted functional repertoire as compared to healthy controls. This was demonstrated by the lower number of distinct avalanche configurations. We focused on the flexibility of the brain activity, since efficient reconfiguration of areas activating above threshold avoids stereotyped, repetitive dynamics, and has been shown to be important for brain functioning 31 . The qualitative restriction in the number of patterns of activation explored by PD patients reflects the effect of pathophysiological changes on the large-scale brain dynamics. A recent EEG study that addressed large-scale dynamics found that the duration, rather than the specific order, of functional states related with symptomatology in patients with Lewy body dementia 32 , which is a disease that shows clinical and pathological similarities to PD 33 . However, the same relationship could not be found in PD patients 32 . Our results add to these findings, indicating that the functional repertoire itself is relevant in PD and related to the clinical picture. In our study, we also tested if the restricted functional repertoire was a result of a slower rate of switching (i.e. from above to below threshold, or vice versa) of each region. Interestingly, we found a higher switching rate in PD compared to the healthy controls, suggesting that the restriction of the functional repertoire in PD is qualitative in nature, and not due to a slowing of the regional rate of threshold crossing. The length of the analyzed time series did not differ significantly between the two groups, thus this result could not be explained by the duration of the acquisitions. Importantly, our results replicated when frequency-specific signals were analysed, showing that the slowing of activity that is normally seen in PD, such as the typical increase in theta and low alpha power 34 , cannot explain the observed differences in the size of the functional repertoire. In fact, even if these frequency bands are filtered out, the functional repertoire of patients is restricted.
We analyzed shared and group-specific avalanche patterns. We classified any pattern as 'shared' if it occurred in both groups, and as 'group-specific' if it occurred in either group, in order to test whether some regions are specifically important in defining pathological activity patterns. The basal ganglia belonged to group-specific patterns more often than expected by chance, suggesting that the network of regions that is dynamically connected with the basal ganglia is altered in Parkinson's disease. This finding highlights the importance of a network www.nature.com/scientificreports/ perspective on brain activity, whereby the basal ganglia show their role in defining the whole-brain dynamics and functional repertoire. This evidence speaks to a vast literature on the basal ganglia as part of a distributed forebrain network with a prominent role in innate behavioral routines as well as learning (retrieving) new (previously acquired) functional repertoires 35 . Previous evidence has shown that extreme specificity is necessary in the connections between the basal ganglia, thalamus and cortex 36 , and that the basal ganglia play a prominent role in the fine tuning of such connections 37 . Our results suggest that the functional role of the basal ganglia might contribute to appropriate dynamical recruitment of brain regions. Importantly, we tested if the restriction of the functional repertoire related to behavioral outcomes. In fact, according to our hypothesis, we expect that impaired brain dynamics at rest would underpin reduced capability to accomplish a variety of complex tasks. The Unified Parkinson's disease rating scale-III (UPDRS-III) is widely used to assess PD patients, and it provides a quick account of motor symptoms severity in PD. In line with our hypothesis, an inverse relationship was evident between the restriction of the functional repertoire and the UPDRS-III. This evidence supports the interpretation of the restricted functional repertoire as pathological, and points towards the idea that an impaired regulation of the brain dynamics by the basal ganglia underlies clinical impairment in PD. In the same line of thinking, the positive correlations that were present between performance on specific cognitive domains and the size of the functional repertoire, confirms that flexible activity is important for performance on a number of behavioral tasks, not limited to the motor domain. Our findings might also be related to the idea of "dedifferentiation" in aging. Aging is a process that causes less efficient dynamics in the brain. A number of studies showed that the same task prompts localized activations in the young, and much broader activations in healthy elderly 38 . This evidence has been interpreted in terms of a compensatory mechanisms, whereby elderly subjects need to recruit more regions to maintain the functional outcomes. Our results would be in line with an absence of such a compensatory mechanism in PD patients. However, future studies involving different tasks are required to elucidate this potential relationship between the size of functional repertoire and dedifferentiation. Finally, our results also highlight the importance of future longitudinal studies to test the performance of the size of the functional repertoire as a potential biomarker of disease progression.
With regard to the mechanisms underpinning the lack of flexibility, these might involve perturbations in the ability to synchronize different areas 27 . Since PD has often been associated with excessive synchronization 39 , we expected to find a reduction in the size of functional repertoire for patients with a more pronounced hypersynchronization in the beta band, which is classically involved in PD 30 . Our results confirmed that PD patients display hypersynchronization, as well as the predicted negative relation between the number of different patterns of activation and the overall synchronization in the beta band.
Electrophysiological studies in animal models of PD 29 and in patients with PD have shown that dopaminergic depletion can increase oscillatory activity in the thalamo-cortical circuitry. Furthermore, while this activity has been described in multiple frequency bands, the beta band is specifically relevant to PD 30 . The extent of such aberrant beta synchronization has been related to clinical disability, and L-DOPA has been shown to partly revert this phenomenon and relieve symptoms 30 . Accordingly, EEG data from PD patients showed that the lack of dopamine increases coupling strength in the basal ganglia circuitry 39 . Within this framework, our results indicate that a hypersynchronized regime is deleterious as it reduces the flexibility of brain activity. In fact, excessive synchronization reduces the variability of the behaviour of the system, hence narrowing the number of functional states that are readily accessible 27,40 . Hypersynchronization and hyperconnected topology have been associated with a variety of neurological diseases 41,42 , as well as in preclinical conditions that carry an increased risk of neurodegeneration 43 , hence investigations of the functional repertoires in these diseases and conditions could help characterizing these conditions. Some limitations of our work should be considered when interpreting our results. First, we utilised the automated anatomicl labelling (AAL) atlas when reconstructing the neuronal activity, which is well suited to MEG studies but has a somewhat coarse spatial resolution. Taking into account that this work is based on MEG signals, it is important to consider that more fine-grained atlases might go over the nominal resolution of the source-reconstructed signal, which could lead to spurious results. Interestingly, parcellations that are based on MEG resting-state data contain roughly the same number of parcels as the AAL atlas 44 . Hence, we believe that this choice of atlas yields a reasonable compromise between obtaining high spatio-temporal resolution while avoiding an overestimation of the coactivations that might be induced by atlas with a too high spatial resolution. Furthermore, we chose to focus our main analysis on 90 regions, excluding the cerebellum, as the parcellation of the cerebellum in the AAL atlas is particularly fine grained, the posterior fossa is prone to artefact, and MEG does not provide high spatial resolution. However, the role of the cerebellum is undoubtedly important in PD 45 . For this reason, we also repeated the analysis with all 116 AAL regions, including the cerebellum. All the significant group differences were confirmed, proving the robustness of the findings. Furthermore, we chose to include the basal ganglia in the analysis, given the paramount importance of their role in PD. Nonetheless, it is worth to consider that the sensitivity of MEG, and thereby the quality of the source reconstruction, tends to decay with the depth of the analysed brain structure. Hence, the analysis involving the basal ganglia is based on a reconstructed signal that may have suboptimal resolution (as compared to the cortical data). However, it has recently been shown that magnetoencephalography can indeed detect the activity of deep brain sources 46,47 .
In conclusion, our results show that the brain of PD patients is less flexible, and a more stereotyped brain activity impairs the ability of the brain to perform complex tasks, and is related to hypersynchronization. Crucially, the size of the functional repertoire is proportional to the observed clinical disability. Within this framework, the main symptoms of PD patients can be explained, and the known role the basal ganglia play in defining PD pathophysiology naturally emerges in the marked involvement in pathological configurations. Our findings also indicate that the mechanisms underlying impaired flexibility of brain activity relates with the observed clinical phenotype observed in PD, and to hypersynchronization. This also provides a new interpretation of PD pathophysiology in which beta hypersynchronization is associated with reduced flexibility of brain dynamics.  Italy). Inclusion criteria were: (a) PD onset after the age of 40 years, to exclude early onset parkinsonism; (b) a modified Hoehn and Yahr (H&Y) stage ≤ 2.5. Exclusion criteria were: (a) dementia associated with PD according to consensus criteria 49 ; (b) relevant cognitive impairment, as evidenced by age-and education-adjusted MoCA score lower than or equal to the Italian cut-off score 50 ; (c) any other neurological disorder or clinically significant or unstable medical condition (Table 1). Disease severity was assessed using the Hoehn and Yahr (H&Y) stages 51 and the UPDRS III 52 . Motor clinical assessment was performed in "off-state" (off-medication overnight). Levodopa equivalent daily dose (LEDD) was calculated for both dopamine agonists (LEDD-DA) and dopamine agonists + L-dopa (total LEDD) 53 . Global cognition was assessed by means of Montreal Cognitive Assessment (MoCA) 54 . MoCA consists of 12 subtasks exploring the following cognitive domains: (1) memory (score range 0-5), assessed by means of delayed recall of five nouns, after two verbal presentations; (2) visuospatial abilities (score range 0-4), assessed by a clock-drawing task (3 points) and by copying of a cube (1 point); (3) executive functions (score range 0-4), assessed by means of a brief version of the Trail Making B task (1 point). All patients were right-handed. All participants signed informed consent. The study was approved by the Local Ethics Committee of University of Naples "L. Vanvitelli" and was conducted in accordance to the Declaration of Helsinki (Table 1).

MEG acquisition.
Magnetoencephalographic data were acquired in a 163-magnetometers MEG system 55 , with 9 reference sensors, located in a magnetically shielded room (AtB Biomag UG, Ulm, Germany). The position of four position coils and four reference points (nasion, right and left pre-auricular and apex) were digitized before acquisition, using Fastrak (Polhemus). MEG data were acquired during two eyes closed resting-state segments each 2.5 min long. Participants were in off-state, and were requested to relax with eyes closed and not to think of anything in particular. Instructions were delivered immediately prior to each recording via an intercom. Head position was recorded at the start of each recording segment. After an anti-aliasing filter, the data were sampled at 1024 Hz. A 4th-order Butterworth IIR band-pass filter was then applied to remove components below 0.5 and above 48 Hz. The filter was implemented offline using MatLab scripts within the Fieldtrip toolbox 2014 56 . Electrocardiogram (ECG) and electrooculogram (EOG) data were also recorded.
MRI acquisition. MR images were acquired on a 3-T scanner equipped with an 8-channel parallel head coil (General Electric Healthcare, Milwaukee, WI, USA) either after, or a minimum of 21 days (but not more than one month) before, the MEG recording. Three-dimensional T1-weighted images (gradient-echo sequence Inversion Recovery prepared Fast Spoiled Gradient Recalled-echo, time repetition = 6988 ms, TI = 1100 ms, TE = 3.9 ms, flip angle = 10, voxel size = 1 × 1 × 1.2 mm 3 ) were acquired.

Preprocessing.
A principal component analysis (PCA) was performed to reduce the environmental magnetic noise 57 . Specifically, the filter was obtained by orthogonalizing the reference signals to obtain a base, projecting the signals from the brain sensors on this noise-base, and subsequently removing these projections www.nature.com/scientificreports/ in order to obtain clean data 58 . We adopted the PCA filtering implementation available within the Fieldtrip Toolbox 56 . The noisy segments of acquisition were identified through visual inspection of the whole dataset by an experienced rater (RR). On average, 130 ± 2 channels were used. After that, Independent component analysis (ICA) 59 was also performed to eliminate ECG (typically 1-2 two components) and EOG (0-1 components) contributions to the MEG signals.
Source reconstruction. Source reconstruction of channel data was performed using a beamforming procedure implemented in the Fieldtrip toolbox 56 . First, the subject's fiducial points were used to co-register the MEG data to the native subject-specific MRI. Second, using a single shell volume conduction model 60 and an equivalent current dipole source model, a Linearly Constrained Minimum Variance (LCMV) beamformer 61 , based on the whole pre-processed, broad-band data, was used to reconstruct time series related to the centroids of 116 regions-of-interest (ROIs), derived from the Automated Anatomical Labeling (AAL) atlas 62,63 . Both the atlas and the MRI were aligned to the head coordinates. We focussed on the first 90 ROIs, excluding those in the cerebellum given the low reliability of the reconstructed signal in this region. For each source, we projected the time series along the dipole direction that explained most variance by means of singular value decomposition (SVD). For each subject, we visually inspected the source-space data to check that no remaining artefact was present at the source-level.

Analysis of dynamics. Neuronal avalanche and avalanche configuration.
To quantify spatio-temporal fluctuations of activity, we first estimated "Neuronal Avalanches". To start, each of the 90 source-reconstructed signals were z-transformed (Fig. 1A). Subsequently, each time series was thresholded according to a cut-off of 3 standard deviations (i.e. z = 3). However, the results are not strictly dependent upon the choice of this threshold.
In fact, we repeated the analyses setting the threshold to z = 2.5 and z = 3.5, and obtained similar results (see Supplementary Fig. 2). We defined a neuronal avalanche as an event in which large fluctuations of activity are present, starting when at least one region becomes active (i.e., above the threshold of (|z|> 3) and continuing as long as any region remains above threshold (Fig. 1B). Note that the position of the regions in space are not taken into account, and hence two coactivated areas (in time) need not be adjacent and may even be located in different hemispheres. These analyses require the time series to be binned. To select a suitable bin length, we computed the branching ratio 2, 64 , σ, as follows: for each time bin duration, for each subject, for each avalanche, the (geometrically) averaged ratio of the number of events (activations) between the subsequent time bin and that in the current time bin was calculated as, where σ i is the branching parameter of the i-th avalanche in the dataset, N bin is the total number of bins in the i-th avalanche, n events j is the total number of events in the j-th bin. We then (geometrically) averaged the results over all avalanches 65 , where N aval is the total number of avalanche in each participant's dataset. In branching processes, a branching ratio of σ = 1 indicates critical processes with activity that is highly variable and nearly sustained, σ < 1 indicates subcritical processes in which the activity quickly dies out, and σ > 1 indicates supercritical processes in which the activity increases. The bin length equal to three samples yielded a critical process with σ = 1, justifying the use of the term "avalanche" for these events. This means that each bin is obtained from three time-points of the binarized time-series. However, the robustness of the results to changes in this exact bin length was also investigated, and the main results were not affected (Supplementary Information). That is, different binnings yielded branching ratios very close to one, and the statistical differences in the exploration of the functional repertoire were similar for each binning (see Supplementary Fig. 1).For each avalanche, an "avalanche configuration" was defined as the set of all areas that were above threshold at some point during the avalanche (see Fig. 1C). Note that the definition of avalanche configuration implies a loss of information in terms of the temporal structure within each avalanche.
Functional repertoire. For each subject, we estimated the functional repertoire, defined as the number of unique avalanche configurations that was expressed during the recording. "Unique" avalanche configurations that each avalanche pattern only counts once towards the size of the functional repertoire (i.e. it does not matter if a given avalanche configuration appear only once or multiple times, as only the number of different configurations contribute to the functional repertoire). Before comparing the functional repertoires between groups, the duration of the acquisitions was compared between the two groups, as the acquisition duration could affect the estimate of the functional repertoire. Furthermore, the analysis was also repeated with the (exact) same amount of data for each participant (67.01 s). To do so, we selected all participant who had at least a 1 min long acquisition (two patients had less and had to be excluded, leading to 38 controls and 37 patients). We randomly selected data segments from the whole recordings. The results from analysis of these data confirmed our initial analysis ( Supplementary Fig. 4).
(1) www.nature.com/scientificreports/ Switching between states. We defined a "switch" as the occurrence of a crossing of the threshold in either direction) between two consecutive time-bins. The switch rate (number of switches over duration), averaged over areas, was computed for each individual.

Analysis of similarity.
To estimate and compare the variability of avalanche configurations between groups, we first computed the similarities between each configuration, yielding one distance matrix per group 2 . In each matrix, rows and columns are equal to unique avalanche patterns, while the entries are the Hamming distances (the number of regions that differ-i.e. above or below threshold-in two given patterns). Statistical comparison of the similarity matrices belonging to the two groups was achieved by permutation testing (see below for details).
Quantifying the influence of specific regions on avalanche configurations. We split the total functional repertoire in two groups: configurations that occurred in both the clinical and control participants ("shared repertoire"), and configurations that were unique to either group ("group-specific repertoire"). We then computed how often each region occurred in both shared and group-specific configurations. The Kolmogorov-Smirnov test was used to compare the two distributions, and resampling then allowed identification of those regions that appeared in group-specific repertoires more often than by chance.

Synchrony estimation.
To estimate synchronization in the beta band, the broadband source-reconstructed data were band-pass filtered with a fourth order Butterworth filter between 13 and 30 Hz. A metric recently developed by our group 66 , the phase linearity measurement (PLM), illustrated in Fig. 5, was then employed to estimate synchronization. The PLM is defined as: where �φ(t) is the interferometric phase (i.e. the phase of the signal resulting from the multiplication of one signal with the complex conjugate of another signal). The PLM was computed using the implementation available in the Fieldtrip toolbox in Matlab. In short, this metric quantifies the synchronization between two time-series using the central peak in the frequency spectrum of the interferometric signal. The more peaked the spectrum of the interferometric signal, the more the two originating signals will be synchronized. The PLM ranges between 0 and 1, is insensitive to volume conduction and grows monotonically with synchronization 66 .

Statistical analysis.
To compare age and sex between the two groups we used T-test and Chi-square test, respectively. Permutation testing or Kolmogorov-Smirnov test were performed to compare patients and controls, as appropriate. For permutation testing, the data where permuted 10,000 times, and at each iteration the absolute value of the difference between the two groups was observed, building a null distribution of absolute differences. Finally, the empirical, observed difference was rank-ordered against this distribution, yielding a significance value. All statistical analyses were performed using custom scripts written in Matlab 2018a.

Data availability
The MEG data and the reconstructed avalanches are available upon reasonable request to the corresponding author, conditional on appropriate ethics approval at the local site. Figure 5. Phase linearity measurement, used to estimate synchronization. The time series of each pair of regions have been processed to obtain the interferometric signal. The power spectrum of the interferometric signal was estimated, and the central area underneath the spectrum (depicted in light brown) was compared to the total area. The higher such ratio, the more synchronized the two signals.