Neuronal Avalanches to Study the Coordination of Large-Scale Brain Activity: Application to Rett Syndrome

Many complex systems, such as the brain, display large-scale coordinated interactions that create ordered patterns. Classically, such patterns have been studied using the framework of criticality, i.e., at a transition point between two qualitatively distinct patterns. This kind of system is generally characterized by a scale-invariant organization, in space and time, optimally described by a power-law distribution whose slope is quantified by an exponent α. The dynamics of these systems is characterized by alternating periods of activations, called avalanches, with quiescent periods. To maximize its efficiency, the system must find a trade-off between its stability and ease of propagation of activation, which is achieved by a branching process. It is quantified by a branching parameter σ defined as the average ratio between the number of activations in consecutive time bins. The brain is itself a complex system and its activity can be described as a series of neuronal avalanches. It is known that critical aspects of brain dynamics are modeled with a branching parameter σ = , and the neuronal avalanches distribution fits well with a power law distribution exponent α = -3/2. The aim of our work was to study a self-organized criticality system in which there was a change in neuronal circuits due to genetic causes. To this end, we have compared the characteristics of neuronal avalanches in a group of 10 patients affected by Rett syndrome, during an open-eye resting-state condition estimated using magnetoencephalography, with respect to 10 healthy subjects. The analysis was performed both in broadband and in the five canonical frequency bands. We found, for both groups, a branching parameter close to 1. In this critical condition, Rett patients show a lower distribution parameter α in the delta and broadband. These results suggest that the large-scale coordination of activity occurs to a lesser extent in RTT patients.


INTRODUCTION
Rett syndrome (RTT) is a severe neurodevelopmental disorder, characterized by clinically normal development for the first 12-18 months of life, when an overall arrest and regression of the psychomotor development begins. The patients lose the acquired verbal skills, the normal motor hand function is replaced with stereotyped movements, and often ataxia is present, which seriously compromises motor skills. Intellectual disability, epilepsy, autonomic dysfunction, breathing abnormalities, anxiety, and orthopedic problems are the other most important clinical findings in RTT (Neul et al., 2010). In 1999, Amir et al. (1999) identified in a mutation in the X-linked MECP2 gene, encoding the methyl CpG-binding protein 2 (MeCP2), the most common cause of RTT.
The dysfunction of MeCP2 provokes multiple effects, such as impaired neuronal maturation, altered GABAergic signaling, and, more importantly, a local imbalance between neuronal excitation and inhibition at the circuit level (Pohodich and Zoghbi, 2015). Importantly, all these molecular and circuital evidence are scarcely linked to the clinically evident deficits in higher cognitive functions. In fact, effective coordination of largescale activity is considered necessary for the emergence of highlevel cognitive abilities. All these features make RTT a potential model for studying how an alteration of the normal development process interferes with the activity of brain networks.
Consequently, we borrowed from physics some tools that allow the description of large-scale activity. It has been shown that the healthy brain has a higher number of configurations compared to the unhealthy brain . Importantly, it has been argued that flexibility is possible given that the brain sets itself to operate near a critical state (Beggs and Plenz, 2003). Bak et al. (1987) first described the idea of "selforganized criticality" (SOC) defined as a system autonomously evolving to a critical state regardless of the initial condition. As he defines, "the critical state is an attractor for the dynamics of a system." In the critical state, the activity can spread across the system in a controlled way, avoiding that activations either spread uncontrollably or quickly die out (Zapperi et al., 1995; Abbreviations: HS, healthy subjects; MECP2, methyl CpG-binding protein 2; MEG, magnetoencephalography; ROIs, regions of interests; RTT, Rett syndrome; SOC, self-organized criticality. Harris, 2002). Such state is classically quantified by a branching parameter σ defined as the average ratio between the number of activations in consecutive time bins (i.e., the sample units) (De Carvalho and Prado, 2000). More specifically, when σ < 1, the system is in a subcritical condition in which the activations quickly fade out. On the other hand, when σ > 1, the opposite happens, and the system is in a supercritical state whereby the activations spread across the entire system, leading to unstable runaway activations. Finally, when σ = 1, the system is in a critical state with balanced, fine-tuned propagation of the activity (Bak et al., 1987;Harris, 2002).
SOC systems are generally characterized by a scale-invariant organization in space and time. In other words, the fluctuations that occur are similar in all scales of time and space, and they can be optimally described by fractal metrics, typically a power law distribution whose slope quantifies the amplitude distribution of such fluctuations (Paczuski et al., 1996;Newman, 2007;Clauset et al., 2009).
Importantly, power law distributions do not prove criticality (Clauset et al., 2009). However, they do provide a statistical measure of the occurrence of large-scale coordinated activity, perhaps providing a useful tool to probe hypothesis about the nature of high-level cognitive deficits, such as those observed in Rett syndrome.
Generally speaking, SOC can be found in very different types of phenomena, like in many natural systems such as earthquakes (Gutenberg and Richter, 1956) or forest fires (Malamud et al., 1998), in which when a unit exceeds a threshold, other units follow it, thus giving rise to a cascade (or avalanche) that spreads across the entire system. The presence of power law-distributed bursts of activations, called avalanches, alternating with quiescent periods, characterizes the dynamics of SOC systems (Lombardi et al., 2016) likely maximizing the efficiency of communication among its elements. However, it is to be considered that the analysis of avalanches, with a power law distribution, varies in heterogeneous/disordered systems and that to get information on the non-equilibrium dynamics, it would be useful to also analyze the crackling noise (Salje and Dahmen, 2014).
In this article, we hypothesize that the genetic defect causing RTT impairs the whole-brain dynamic and relates to clinical symptoms. To test our hypothesis, we compared the distributions of the size of the neuronal avalanches, in order to capture the critical features of brain activity, in 10 patients with RTT and  10n matched healthy subjects (HS). We expect that large-scale activation would be occurring less in patients showing marked deficits in higher cognitive functions. We based our analyses on source-reconstructed magnetoencephalographic (MEG) data acquired during resting state.

Participants and Clinical Assessment
Participants were recruited from the Department of Translational Medical Sciences, Child Neuropsychiatry in "Federico II" University of Naples, Italy. We studied 10 female patients with clinical diagnosis of RTT (age 24.30 ± 8 years) based on the Neul revised criteria (Neul et al., 2010) and confirmed by mutation in the MECP2 gene, and 10 healthy female individuals (HS) (age 26.10 ± 6.84 years). This study complied with the Declaration of Helsinki and was approved by the local ethics committee. Written informed consent has been granted by all participants (or their legal guardians).

Acquisition
The data were acquired using a MEG system equipped by 163 magnetometers SQUID (Superconducting Quantum Interference Device) , placed in a magnetically shielded room (AtB Biomag, Ulm, Germany), in order to  For brevity, we reported only delta, alpha, and broadbands.
FIGURE 2 | Size distributions of avalanches with the relative power law fitted, for each RTT patient (first line) and for each HS (second line). The colored lines represent the size distributions of avalanches for each subject (RTT patients in the top row, Healthy Subjects in the bottom row). The black bold lines represent the fitted power laws. For brevity, we reported only Delta ( t = 5), Alpha ( t = 1) and Broad ( t = 2) bands. Specifically, the Delta band and the Broad band reached statistical significance, whereas the Alpha band was reported for comparison, as an example, as it does not reach significance (such as Theta, Beta and Gamma). The reported results refer to the threshold SD = 3.
reduce the external noise. Of all squids, 154 are positioned to be as close as possible to the head of the subject, while the remaining ones are more distant so as to measure environmental noise (reference magnetometers). All the subjects of each group underwent a 7-min resting-state MEG acquisition with open eyes, divided in two segments. To evaluate the right  position of the head under the helmet, we used Fastrak (Polhemus R ) to acquire the position of four coils (attached to the head) and of four anatomical landmarks (nasion, right and left pre-auricular points, and vertex of the head). Before each acquisition segment, the head position in the helmet was obtained. During the acquisition, two electrodes for the electrocardiogram and two for the electrooculogram (Gross et al., 2013) were also acquired.

Preprocessing
As previously described Rucco et al., 2019), the MEG signals, after an anti-aliasing filter, were acquired with a sampling frequency of 1,024 Hz. A fourth-order Butterworth IIR band-pass filter in the 0.5-to 100-Hz band was subsequently applied to the acquired signals. Environmental noise, measured by reference magnetometers, was removed by using the principal component analysis. MEG data were cleaned of physiological artifacts, such as eye blinking and heart activity, by means of independent component analysis (Sorriso et al., 2019). Visual inspection was used for identification of noisy channels. For all the preprocessing steps, we used the FieldTrip toolbox (Oostenveld et al., 2011).

Source Reconstruction
Time  corresponding to the cerebellum because of their low reliability in MEG . Hence, we considered a total of 90 ROIs. We have resampled the source-space time series at 512 Hz.

Signal Discretization
For each ROI, the time series were discretized with five different time bin durations t, each one multiple of t min (19 ms), corresponding to the sampling period. An event was identified by a positive or a negative excursion of an area, in a bin, beyond a threshold, defined as ±3 SD of the signal amplitude, which was a tradeoff between a lower threshold (leading to the detection of more spurious noise events in addition to real ones) and a higher threshold (missing real events). Indeed, MEG systems rely on squid sensors, which measure the magnetic field generated by brain activity. However, such sensors cannot identify the field direction. For this reason, both positive and negative excursions were picked up without distinction. As a result, the event was identified when the absolute value of this excursion exceeded the chosen threshold.

Branching Parameter
As previously described (Beggs and Plenz, 2003;Shriki et al., 2013), an avalanche is defined as a sequence of contiguous time bins starting when at least one ROI is active and ending when all regions are inactive. The number of events in all ROIs in an avalanche corresponds to its size.
For each subject and for each time bin size, the branching parameter σ was estimated by calculating, for each avalanche, the averaged (over all the time bins) ratio of the number of events between the subsequent time bin (descendants) and that in the current time bin (ancestors), and then averaging it over all the cascades (Bak et al., 1987). More specifically: where σ i is the branching parameter of the i-th avalanche in the dataset, N bin is the total amount of bins in the i-th avalanche, n events (j) is the total number of events active in the j-th bin, and N aval is the total number of avalanche in the dataset.

Neuronal Avalanches and Power Law Parameter
The statistical distribution of the avalanche size x was assumed to be described by a power law function: where C α is a normalization factor. The parameters x min was set to 1 (the minimal avalanche size), and x max was set to 1.5 times the total number of ROIs included in the analysis. We used truncated power law to take finite size effect into account. The distribution parameter α was estimated from the data by implementing a maximum likelihood algorithm (Kay, 1993). When the branching parameter σ is equal to 1, the system is in a critical state. In other words, there is a long distance propagation of neuronal activity without runaway excitation. For this reason, we calculated at which time bin, for each frequency band, the branching parameter (as average distance) was the closest to 1 for both groups. Such choice was made in order to FIGURE 8 | Duration versus pause for both groups for delta ( t = 5), alpha ( t = 1), and broad ( t = 2) bands varying the threshold (1 to 3 SD). have the two groups in the same condition, in order to allow a meaningful comparison of the α parameter. Furthermore, to evaluate the robustness of our data, we split our dataset both in two and in three segments, and for each of them, we calculated the sigma parameter and his variance. Figure 1 shows the data analysis pipeline.

Goodness of Fit
We performed the goodness of fit (GOF) via the evaluation of Hellinger distance (HD) (Oosterhoff and van Zwet, 2012) between the distribution of real data and the distribution of fitted data, in order to quantify the similarity between two probability distributions, as follows.
where P is the probability mass function of real data, and Q is the probability mass function of fitted data. The HD is in between 0 and 1; and the lower the HD, the better the fit to real data (ideally, HD = 0). Furthermore, to assess the GOF overall RTT patients and HS, we provided the empirical cumulative distribution function (ECDF) of HD.

Statistical Analysis
We compared the distribution parameter α in patients and controls using permutation testing (Nichols and Holmes, 2002), for each frequency band. Specifically, at each iteration, each subject is randomly assigned to one of the two groups and then the difference between the averages of the two groups was computed. The group assignment is permuted 104 times, obtaining the null distribution of group differences, which was used to define the statistical significance of the observed difference between patients and controls. All the analyses were performed at a significance level of 0.05, using Matlab R2017a (MathWorks R ) environment.

RESULTS
As described in the Methods section, we discretized the ROI signals as a train of events, each of them representing the signal exceeding the threshold of ±3 SD.
First, we checked the critical condition in both groups, in all frequency bands, for all time bins.
We found, for both groups, a branching parameter σ ≈ 1, in delta band for bin t = 5, in theta and alpha band for bin t = 1, in beta-, gamma-, and broadband for bin t = 2, as reported in Table 1. Moreover, to demonstrate the robustness of the estimation of the branching parameter, we estimated the variance around sigma both when we split our dataset into two segments ( Table 2) and when we divided it into three segments (Table 3).
To demonstrate the strength of our data, we reported, for the chosen threshold SD = 3, the size distribution fitted with a power law (Figure 2), duration distributions (Figure 3), and average size versus duration distributions (Figure 4) for each subject, in all time bins, for both groups separately (for brevity, we reported only delta, alpha, and broadbands). Again, we provided the effect of different thresholds (1 to 3 SD) on our data in size distributions (Figure 5), duration distribution (Figure 6), and average size versus duration distributions (Figure 7).
Furthermore, to test that our results might not be biased by spatial subsampling, we demonstrate the presence of scale time separation (Levina and Priesemann, 2017), as reported in Figure 8.
Last, to quantify the robustness of power law fitting, we evaluated the ECDF of HD for both groups, for all time bins in delta, alpha, and broadbands, and as reported in Figure 9, the highest HD is less than 0.3.
Finally, a comparison of the α parameter between both groups was carried out, highlighting significant differences in the delta band (p = 0.0423) and in broadband (p = 0.0174), as shown in Figure 10. In particular, RTT patients invariably showed a lower distribution parameter α in both frequency bands, for all time bins. No significant differences were found in other frequency bands.

DISCUSSION
In this work, we set out to test the hypothesis that genetic mutations related to RTT would induce a reorganization of the whole-brain activity.
In particular, we characterized the distribution of neuronal avalanches in a group of 10 patients with Rett syndrome and 10 matched controls, using MEG, in an open-eye resting-state condition. The analysis was performed both in broadband and in the five canonical frequency bands. The study originated from recent evidence showing that the dynamics of resting-state brain activity, measured using MEG, produces scale-invariant neuronal avalanches, suggesting that the critical state is a physiological condition that is (presumably) optimal to the brain functioning. Consequently, the deviations from this optimally tuned configuration might convey the effects of the disease on large-scale coordinated activity (Meisel et al., 2012).
Our data are in line with recent findings showing that, in humans, the dynamics of resting-state brain activity is well described by a branching process where the sigma parameter is close to 1 (Meisel et al., 2013;Shriki et al., 2013).
Besides these confirmatory results, the main goal of this work is to study the frequency-specific differences in functional dynamics in RTT patients with respect to a control group. Our data showed that, around the critical state, the distribution size of neuronal avalanches, for both groups, obeys a power law. Interestingly, RTT patients show a lower value of the distribution parameter α both in delta and broadband, compared to controls. Furthermore, the exponents of the power law distributions vary greatly per each frequency band. This result is in line with Thompson et al., which suggests that resting-state connectivity is a frequency-dependent phenomenon (Thompson and Fransson, 2015).
The aforementioned results collectively suggest a global rearrangement in the brain functional dynamics of Rett patients. In fact, big avalanches capture the presence of widespread coordinated activations, with many different brain areas activating in a complex, meaningful sequence. The lack of large avalanches in patients might be a sign of the inability of the brain to coordinate a sufficient number of areas in sequence. One might speculate that the lack of large patterns of structured activity might be capturing the effects of the circuital changes induced by the genetic mutations present in RTT. This hypothesis might be compatible with the work by Armstrong et al., which claims that whole-brain structural changes are related to small, low-density neurons that have less dendritic branching and spine density (Armstrong et al., 1995). Furthermore, in the RTT syndrome, the dendrites of pyramidal neurons in the motor and frontal cortices are considerably shorter than in HS (Armstrong et al., 1995).
While our results are not conclusive about the brain being operating at criticality, they appear to be compatible with such hypothesis. This might be relevant given that it is known that critical dynamics facilitates optimization of processes such as information transfer and storage capacity. In particular, SOC involves a number advantages in information processing: (i) optimal information transmission because there is a perfect balance between the propagation of the signal over long distances and the resistance to saturation (Beggs and Plenz, 2003;Bertschinger and Natschläger, 2004;Kinouchi and Copelli, 2006); (ii) maximum mnemonic repertoire size, as the number of significant repetitive avalanche patterns is maximized (Haldeman and Beggs, 2005;Shew et al., 2011;Fagerholm et al., 2016); and (iii) maximized computational performance, as near the critical point, networks have better performances in different computational tasks, compared to networks with subcritical or supercritical dynamics (Bertschinger and Natschläger, 2004).

CONCLUSION
Concluding, in our study we analyzed brain dynamics during eye-open resting state, measuring the size distribution of neuronal avalanches. We found a lower distribution parameter α in RTT patients in the delta frequency band and in the broadband, compared to controls. Such finding means that patients display a larger number of small avalanches and a lower number of big avalanches. These results can be explained in terms of reduced long-range coordination of neuronal activity across the brain in the pathological group. Such altered brain dynamics is likely to be suboptimal, and might underpin some of the symptomatology observed in Rett syndrome. Additionally, our results suggest that neuronal avalanches could be an innovative and advanced method to study the dynamics of brain activity. In fact, despite the rising interest in what constitutes the normal cortical dynamics in healthy humans, the nature of the alterations induced by pathological processes remains a major question in neuroscience. Interestingly, the observation that cortical networks optimize information processing when they are in a critical state might be exploited to identify subtle, preclinical brain dysfunction. Furthermore, the study of the critical dynamics of the human brain could be enhanced by comparing normal resting avalanches with evoked states or pharmacologically modified states and assessing the sensitivity (or performing a profile analysis) of the system with respect to the branching parameter values. Finally, such approaches are inherently multivariate, taking into account the activity of all brain areas at once, hence providing solid theoretical ground to the study of emerging, holistic properties of the brain.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by "COMITATO ETICO UNIVERSITA' FEDERICO II" n • protocol: 392/19. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
RR collected the sample, performed the MEG recordings, preprocessed the MEG data, contributed to MEG data analysis and statistical analysis, wrote the manuscript, and prepared the figures. PB collected the sample, contributed to clinical data analysis, and wrote the manuscript. AL performed the MEG recordings and preprocessed the MEG data. FB collaborated with the MEG data analysis and with the statistical analysis. MP performed the MEG recordings and preprocessed the MEG data. AP performed the MEG recordings and preprocessed the MEG data. CB collected the sample. CG, LM, and GS supervised the study. PS contributed to interpreting the results and critically revised the article. All authors read and approved the final version of the manuscript.

FUNDING
This work has been partly supported by the University of Napoli Parthenope under the 'Bando di ricerca competitiva di ateneo per il triennio 2016-2018 progetto: implementazione di studi di connettività cerebrale mediante magnetoencefalografia in ambito multidisciplinare'.