Introduction

Brain functioning requires efficient reconfiguration of patterns of activity. This flexibility is essential for the coordinated engagement of brain regions, which underlies complex behaviours. Moreover, dynamic patterns of activity reflect the coordinated engagement of different systems within the brain1. Accordingly, a large number of distinct patterns of activity (“functional repertoire”) indicates flexible dynamics. Growing evidence also indicates that the maximum number of spatio-temporal patterns of activity occurs at criticality and deviations from this regime incur in a reduction of the number of patterns observed2. Criticality is a highly variable, adaptive and flexible dynamical regime3. It optimizes the capability of storing information4, the efficient transmission of information across the brain3, the response to internal fluctuations5 and the detection of external stimuli6; for review see Ref.7.

The study of temporally resolved patterns of activation has recently drawn significant interest8. Furthermore, it has been shown that they may underpin cognitive and behavioural functions9. In fact, a number of behavioural functions (both motor and cognitive) require flexible dynamics, and brain pathologies might shift the brain to a dynamical regime where reconfigurations are not efficient, hence impairing behavioural functions that need such flexibility. Notably, accounting for temporal variability of brain activity has improved the diagnostic classification of neurodegenerative diseases, indicating the value of time-resolved brain-network properties10.

Parkinson’s disease (PD) is a disabling neurodegenerative disorder. It is characterized by reduced amplitude of movements, and slowing of cognitive processes, imposing major individual and social burden11. The main histopathological finding in PD is a severe nigrostriatal dopamine depletion12. PD has traditionally been regarded predominantly as a motor disease. However, recent clinical and neuroimaging findings have questioned these views13. For example, the first symptoms to appear are often not motor, and the clinical phenotype clearly goes well beyond the motor impairment, with other domains, such as executive functioning, involved14. Structural Magnetic Resonance Imaging (MRI) has also confirmed that the regions involved in PD are more widespread than previously thought15. Although functional magnetic resonance (fMRI) studies have identified impairment in the corticostriatal network in PD patients, the observed changes in functional connectivity extend into many other brain systems16. Similarly, magnetoencephalography (MEG) studies have also reported widespread alterations in multiple functional connections in PD17.

In particular, the multifaceted nature of large-scale brain interactions in PD is highlighted by contrasting reports in the literature, with studies reporting more-connected networks18, 19, less-connected networks20, 21, and a combination of both22, 23. However, the structure of dynamic fluctuations in brain activity or connectivity between regions was not extensively explored yet.

Recent studies, specifically addressing the dynamical changes in PD, demonstrated that patients dwelled in a hyperconnected dynamical regime more frequently as compared to matched controls24, and showed greater network level integration during the medication-OFF compared to the ON condition25. These findings suggest that impaired flexibility might be deleterious and potentially also related to clinical symptomatology.

One candidate mechanism of impaired flexibility derives from alterations in excitation and synchronizability26, 27. Optimal flexibility requires a moderate amount of synchronizability28. In contrast, hyper-synchronization can restrict the variability of brain configurations (functional states), entraining the dynamics onto a limited number of patterns. Conversely, too little synchronization also limits the variability of brain states with insufficient integration27. This framework may be specifically relevant for PD, where hyperactivation and hypersynchronization (primarily in in the beta band) are linked to clinical disability29. Furthermore, administration of L-DOPA can partly revert aberrant hypersynchronization in the beta band and relieve symptoms30.

Drawing on this theory, we here hypothesized that functional activity in PD would be less flexible than in matched controls. This flexibility is estimated based on the variability of the types of patters of activations as measured from source-reconstructed MEG signals. To quantify such variability, we used the idea of “neuronal avalanches”. In this context, neuronal avalanche refers to a burst of large-scale activations (i.e. a sequence in time of regions that go above threshold). These bursts of activations are not stereotyped, but rather they often change in time. This means that different combinations occur. In this work, we define the functional repertoire as the number of unique combinations (a unique combination means that the spatial pattern of the avalanche does not overlap fully with the spatial pattern of any other avalanche—repetitions are discarded), and use it as a measure of flexibility of the dynamics (see Fig. 1). Finally, we hypothesized that the basal ganglia would play a particularly important role in determining the occurrence of a reduced number of patterns. Furthermore, we also expected greater symptom expression to be associated with a more restricted functional repertoire.

Figure 1
figure 1

Schematic representation of an avalanche pattern. (A) Reconstructed time series (z-scores). The red dashed lines indicate the threshold to define activation. (B) 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. The figure was done using MatlabR2019a and Inkscape 1.0.

Full size image

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–30 Hz) 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).

Figure 2
figure 2

Functional repertoire computed from high spatiotemporal MEG data. (A) The number of unique avalanche configurations in PD and healthy controls (HC). The central mark in the box indicates the median, and the edges of the box the 25th and 75th percentiles. The whiskers extend to 3 standard deviations, and the outliers are plotted individually using the '+' symbol. (B) Avalanche configurations for two representative participants (regions above threshold in yellow, regions below threshold in green). On the x-axis, the 90 AAL regions (cortex and basal ganglia). Along the y axes, one can see each of the unique patterns that constitute the functional repertoire of the individual. (C) The number of switches for the two groups. (D) Boolean similarities between avalanche configurations. Rows and colums show all avalanche configurations that are present in each cohort. The entries are the boolean similarities (i.e. yellow indicates higher similarity). The matrix-entries have been re-ordered using the Louvain modularity algorithm (as implemented in the brain connectivity toolbox), to graphically highlight the presence of clusters. The coloured lines (white and red) show the clusters in patients and controls, respectively. The figure was done using MatlabR2019a and Inkscape 1.0.

Full size image
Figure 3
figure 3

Relationship between flexibility of brain dynamics, clinical disability and global synchronization. (A) Individual whole-brain average PLM in PD and controls (B) Correlation between global synchronization, estimated using the phase linearity measurement (PLM) in the beta band (13–30 Hz), and the size of the functional repertoire. (C) Correlation between motor impairment, estimated using the Unified Parkinson disease rating scale—III, and the size of the functional repertoire.

Full size image

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 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).

Figure 4
figure 4

Mapping regional contribution to functional repertoire. The regions that contributed most often to avalanche patterns that were unique to one group. In the bar plot on the left, the width of the bar corresponds to the group difference in frequency with which a region contributed to the avalanche patterns that were unique to a group. The image was done using MatlabR2019a, including BraiNetViewer v. 1.62—http://www.nitrc.org/projects/bnv/), and Inkscape 1.0.

Full size image

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, gamma 30–48 Hz). 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 flexibility7. 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 functioning31. 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 dementia32, which is a disease that shows clinical and pathological similarities to PD33. However, the same relationship could not be found in PD patients32. 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 power34, 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 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 repertoires35. Previous evidence has shown that extreme specificity is necessary in the connections between the basal ganglia, thalamus and cortex36, and that the basal ganglia play a prominent role in the fine tuning of such connections37. 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 elderly38. 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 areas27. Since PD has often been associated with excessive synchronization39, 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 PD30. 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 PD29 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 PD30. 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 symptoms30. Accordingly, EEG data from PD patients showed that the lack of dopamine increases coupling strength in the basal ganglia circuitry39. 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 accessible27, 40. Hypersynchronization and hyperconnected topology have been associated with a variety of neurological diseases41, 42, as well as in preclinical conditions that carry an increased risk of neurodegeneration43, 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 atlas44. 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 PD45. 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 sources46, 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.

Materials and methods

Cohort description

Consecutive early PD patients, diagnosed according to the diagnostic criteria of the UK Parkinson’s Disease Society Brain Bank Diagnostic Criteria48, were recruited at the Movement Disorders Unit of the First Division of Neurology at the University of Campania “Luigi Vanvitelli” (Naples, 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 criteria49; (b) relevant cognitive impairment, as evidenced by age- and education-adjusted MoCA score lower than or equal to the Italian cut-off score50; (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) stages51 and the UPDRS III52. 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).

Table 1 Demographic and clinical features of PD patients.
Full size table

MEG acquisition

Magnetoencephalographic data were acquired in a 163-magnetometers MEG system55, 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 201456. 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 mm3) were acquired.

Preprocessing

A principal component analysis (PCA) was performed to reduce the environmental magnetic noise57. 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 in order to obtain clean data58. We adopted the PCA filtering implementation available within the Fieldtrip Toolbox56. 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 toolbox56. 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 model60 and an equivalent current dipole source model, a Linearly Constrained Minimum Variance (LCMV) beamformer61, 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) atlas62, 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 ratio2, 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,

$$\sigma_{i} = \frac{1}{{N_{bin} - 1}}\mathop \prod \limits_{j = 1}^{{N_{bin} - 1}} \left( {\frac{{n_{events} \left( {j + 1} \right)}}{{n_{events} \left( j \right)}}} \right)^{{\frac{1}{{N_{bin} - 1}}}} ,$$
(1)

where \(\sigma_{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} \left( j \right)\) is the total number of events in the j-th bin. We then (geometrically) averaged the results over all avalanches65,

$$\sigma = \frac{1}{{N_{aval} }}\mathop \prod \limits_{i = 1}^{{N_{aval} }} \left( {\sigma_{i} } \right)^{{\frac{1}{{N_{aval} }}}} ,$$
(2)

where \(N_{aval}\) is the total number of avalanche in each participant’s dataset. In branching processes, a branching ratio of \(\sigma\) = 1 indicates critical processes with activity that is highly variable and nearly sustained, \(\sigma < 1\) indicates subcritical processes in which the activity quickly dies out, and \(\sigma > 1\) indicates supercritical processes in which the activity increases. The bin length equal to three samples yielded a critical process with \(\sigma\) = 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).

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 group2. 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 group66, the phase linearity measurement (PLM), illustrated in Fig. 5, was then employed to estimate synchronization. The PLM is defined as:

$$PLM = \frac{{\mathop \smallint \nolimits_{ - 1}^{1} \left| {\mathop \smallint \nolimits_{0}^{T} e^{i\Delta \phi \left( t \right)} e^{ - i2\pi ft} dt} \right|^{2} df}}{{\mathop \smallint \nolimits_{ - \infty }^{\infty } \left| {\mathop \smallint \nolimits_{0}^{T} e^{i\Delta \phi \left( t \right)} e^{ - i2\pi ft} dt} \right|^{2} df}},$$

where \(\Delta \phi \left( t \right)\) 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 synchronization66.

Figure 5
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.

Full size image

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.