Topological changes of fast large-scale brain dynamics in Mild Cognitive Impairment predict early memory impairment: a resting-state, source reconstructed, magnetoencephalography study

Functional connectivity has been used as a framework to investigate widespread brain interactions underlying cognitive deficits in Mild Cognitive Impairment (MCI). However, many functional connectivity metrics focus on the average of the periodic activities, disregarding the aperiodic bursts of activity (i.e., the neuronal avalanches) characterising the large-scale dynamic activities of the brain. Here, we apply the recently described Avalanche Transition Matrix framework to source-reconstructed magnetoencephalography signals in a cohort of 32 MCI patients and 32 healthy controls (HC) to describe the spatio-temporal features of neuronal avalanches and explore their topological properties. Our results showed that MCI patients showed a more centralised network (as assessed by higher values of the degree divergence and leaf fraction) as compared to HC. Furthermore, we found that the degree divergence (in the theta band) was predictive of hippocampal memory impairment. These findings highlight the role of the changes of aperiodic bursts in clinical conditions and may contribute to a more thorough phenotypical assessment of patients.


Introduction
Mild cognitive impairment (MCI) is a preclinical, transitional stage between healthy aging and dementia.Its prevalence ranges from 6.7% to 25.2% in people older than 60 years old and increases with age and lower education level.Nowadays, MCI is considered as a time frame during which we should act to delay conversion to dementia (Jongsiriyanyong and Limpawattana, 2018).From a clinical standpoint, MCI patients are classified according to type and number of affected cognitive domains.The first distinction is made on the presence or absence of memory impairment, giving rise to the 2 major subtypes of MCI, namely amnestic (aMCI) and nonamnestic MCI, where the former denotes a memory loss, while the latter refers to a cognitive impairment which spares memory and affects specific cognitive domains such as executive functions, attention, visuospatial ability, or language.In this context, aMCI is regarded as a precursor to Alzheimer's dementia (Buldú et al., 2011;Celone et al., 2006), and the atrophy of the hippocampus, the key region for learning and memory acquisition (Burgess et al., 2002), is considered as a marker of the conversion from MCI to Alzheimer's disease (AD).In this regard, the relationship between the hippocampal volume and both the memory and the functional connectivity (FC) efficiency is well established (Apostolova et al., 2006;Jacini et al., 2018;Stoub et al., 2010;Wang et al., 2006;Yavuz et al., 2007).Indeed, converging lines of evidence has shown that hippocampal atrophy is directly related to a progressive cognitive impairment (especially with respect to the memory retrieval impairment) and FC alterations in both MCI and AD (Bai et al., 2009;Cai et al., 2017).
Of note, these clinical forms may evolve over time, so if more than 1 domain is involved (in both aMCI and nonamnestic MCI), the label of "multiple-domain," as opposed to "single domain," is used (Petersen, 2016).For example, as aMCI progresses, other cognitive functions may become compromised, defining a continuum between amnestic and multiple-domain aMCI.Hence, constraining the cognitive deficits in MCI to the malfunctioning of any specific brain region might be misleading, as they might rather stem from the miss-interactions among multiple brain regions (Liu et al., 2012;Minati et al., 2014;Sorrentino et al., 2021bSorrentino et al., , 2018)).
Network theory has been used as a framework to characterize large-scale interactions, in order to pinpoint specific mechanisms underpinning cognitive deficits.To this end, FC, which typically measures the pair-wise statistical dependencies between regional signals (Friston, 2011), has been applied in health and disease (Fornito et al., 2015).The FC changes reported in MCI were nonhomogeneous even within individual studies, as the strength of the connectivity was increased in some brain regions and decreased in others (Contreras et al., 2017;Jacini et al., 2018;Liu et al., 2012;López et al., 2017;López-Sanz et al., 2017).Moreover, MCI participants display a hypersynchronized network, which then desynchronizes if overt disease sets in (Buldú et al., 2011;Celone et al., 2006).Importantly, the reported FC changes failed to convincingly replicate across studies, and new approaches are needed to provide reliable estimates of large-scale brain interactions.
A source of such variability could arise from a bias in most FC studies, which generally assume a stationary condition in brain activity.These studies, by focusing on "average" connectivity over a time interval, disregard the evolution of the coactivations over time (i.e., brain dynamics) (Zalesky et al., 2014).In fact, it was shown that large-scale brain activity is far from stationary, and instead it is characterized by aperiodic, scale-free bursts of activity (in the context of statistical mechanics often referred to as neuronal avalanches), which are indeed expected, given a nonlinear underlying dynamics (Haldeman and Beggs, 2005;Shriki et al., 2013;Tagliazucchi et al., 2012).Hence, the study of time-averaged connectivity might not be optimal to capture such bursty, nonlinear activities.
The dynamics of a healthy brain has been shown to display high flexibility, since the regions that are recruited change at each subsequent avalanche, generating complex spatio-temporal patterns.Hence, the number of avalanche patterns, that is, the functional repertoire, has been used as a proxy for the flexibility of brain activity.Flexible dynamics is maintained in health and, conversely, is disrupted in disease (Palmigiano et al., 2017;Sorrentino et al., 2021a).In fact, smaller functional repertoires occur in neurological diseases and are predictive of both clinical disability and disease progression (Duma et al., 2023;Polverino et al., 2022;Sorrentino et al., 2021a).
In this study, in analogy with previous findings (Polverino et al., 2022;Sorrentino et al., 2021a), we hypothesize that MCI would result in a more stereotyped brain dynamics, as measured by a smaller functional repertoire.Then, we moved on to an in-depth characterization of the spatio-temporal dynamics of avalanches (Sorrentino et al., 2021c(Sorrentino et al., , 2023a)).In fact, perturbations of local activation do not occur randomly but, rather, they spread preferentially across the white-matter bundles (i.e., the structural connectome) (Sorrentino et al., 2021c(Sorrentino et al., , 2022)).As a consequence, not all regions are recruited equally by the avalanches (as the topology of the structural connectomes will tend to attract the functional dynamics to its eigenmodes) (Tewarie et al., 2022).We hypothesize that the reduction in the flexibility of brain dynamics would be mirrored by a rearrangement of the spreading of the perturbations on the large-scale activity.In particular, we predict that, in MCI, the regions that are structurally more central would be recruited more often as compared to healthy controls (HC) (which achieve a more flexible dynamics).Conversely, more peripheral regions would be recruited even less often in MCI participants.The reasoning behind this hypothesis lies in the fact that, in the healthy brain, the spreading of neuronal avalanches is influenced by the structural connectome, but it is not fully determined by it.In fact, healthy participants display a flexible dynamic that varies over time and "departs" from the structural connectome.In MCI, the lack of flexibility would reduce the number of dynamical reconfigurations, resulting in a stronger influence of the underlying structure over the spatio-temporal evolution of neuronal avalanches.Globally, this would result in changes in a rearrangement of the functional topology toward a more "centralized" network.
To test our hypothesis, we performed magnetoencephalography (MEG) recording on 32 MCI patients and 32 HC.Then, the sourcereconstructed brain signals were analyzed to extract neuronal avalanches, operationally defined as an event starting when large fluctuations of activity are present in at least 1 brain region and ending when all the regions return to their baseline.Hence, we define an avalanche pattern as the set of brain regions recruited in each individual avalanche.In turn, the totality of the unique avalanche patterns (i.e., discarding repetitions) occurring over time is defined as the functional repertoire, and the number of such patterns defines its size, which is used to quantify flexibility.Finally, the topological features of brain dynamics were assessed using the recently developed avalanche transition matrix (ATM).The ATMs convey the spatio-temporal trajectories of neuronal avalanches as they spread across the brain and define the probability that 2 brain regions will move away from their baseline activity (i.e., avalanches are occurring) in 2 consecutive time frames (Sorrentino et al., 2021c).Finally, in order to verify a possible relationship between the topological organization of brain dynamics in MCI and its clinical features, we built a multilinear regression model to test if the topological properties of the ATMs can predict the clinical impairment.

Participants
Thirty-two MCI patients (18 males and 14 females; 21 singledomain aMCI and 11 multiple-domain aMCI; mean age 71.31; SD ± 6.83; mean education10.54;SD ± 4.33) were recruited from the Center of Cognitive and Memory Disorders of the Hermitage Capodimonte Clinic in Naples, Italy.All the patients were righthanded and native Italian speakers.The MCI diagnosis was done according to the National Institute on Ageing-Alzheimer's Association criteria (Albert et al., 2011).Inclusion criteria were (1) the absence of neurological or systemic illness that could affect the cognitive status, (2) no contraindications to MRI or MEG recording, and (3) Fazekas score (for both periventricular white matter and deep white matter scores) ≤2.Thirty-two participants (19 males and 13 females) matched for age (69.9 ± 5.61) and education (12.96 ± 4.56) were enrolled as a control group (HC).The cohort characteristics are summarized in Table 1.Both HC and MCI patients underwent neurological examination, MRI scan (which included the estimation of the hippocampal volumes), and MEG recording.The neuropsychological screening (see Table 2) included also the free and cued selective reminding test (FCSRT), which is highly sensitive in detecting hippocampal amnestic deficits (Auriacombe et al., 2010;Grober and Buschke, 1987).In particular, the FCSRT assesses the ability of encoding and retrieval through semantic cues by exploiting and maximizing the learning effect (Dubois et al., 2014).The first part of the test is defined as "immediate recall," and it aims to induce the semantic encoding while the second part (performed after 20 minutes) is defined as "delayed recall."The sum of these 2 measures (i.e., the total recall) distinguishes with high sensitivity and specificity the older adults with dementia from those without it.For both parts of the tests, the recall was first free and then cued (Auriacombe et al., 2010).
The study protocol was approved by the ''Comitato Etico Campania Centro'' (Prot.n.93C.E./Reg.n.14-17OSS) and all participants provided written informed consent in accordance with the Declaration of Helsinki.

MEG acquisition and preprocessing
Data were acquired using a magnetoencephalography (MEG) system composed of 154 magnetometers, superconductive quantum interference device, and 9 reference sensors.The acquisition took place in a magnetically shielded room (ATB, Biomag, ULM, Germany) to reduce external noise.To define the position of the head under the helmet, we used Fastrack (Polhemus), that digitized the position of 4 anatomical landmarks (nasion, right and left preauricular points, and vertex of the head) and the position of 4 reference coils (attached to the head of the subject).Each subject was recorded twice (3.5 minutes each) with a 1-minute break, in resting state with closed eyes.We also recorded the cardiac activity and the eyes movement to remove physiological artifacts.After applying an antialiasing filter, data were sampled at 1024 Hz.Data preprocessing was performed similarly to Liparoti et al. (2021).Briefly, the MEG data were filtered in the band 0.5-48 Hz through the implementation of a fourth-order Butterworth IIR band-pass filter, using the Fieldtrip toolbox in the MATLAB environment.Moreover, a principal component analysis was carried out to reduce environmental noise.Finally, we performed a supervised independent component analysis to remove physiological artifacts from the electrocardiography (1 component per participant) and the electroculography (no component per participant, rarely one) (Romano et al., 2022;Rucco et al., 2022).
the linearly constrained minimum variance (Van Veen et al., 1997) beamformer algorithm, based on the native MRI of each subject.Indeed, an expert operator, visually inspected the individuals' MRI, and individuated the coordinates of the 4 anatomical fiducials that were recorded on the participants' head before the MEG recording (i.e., nasion, left preauricular, right preauricular, vertex).We then matched the reconstruction of the sources with the structural image of the brain of each participant, using parcels defined based on the automatic anatomical labeling atlas.The time series were filtered in the 5 canonical frequency bands: delta (0.5-4 Hz), theta (4-8 Hz), alpha (8-13 Hz), beta (13-30 Hz), and gamma (30-48 Hz).

Analysis of brain dynamics 2.5.1. Neuronal avalanches and branching parameter
We estimated the "Neuronal Avalanches" to quantify the spatiotemporal fluctuations of brain activity.A neuronal avalanche is defined as an event starting when massive fluctuations of brain activity are present in at least 1 ROI and ending when all the ROIs return to their usual activity (Sorrentino et al., 2021a).Firstly, we calculated the z-score of each ROIs' time series, and then each time series was thresholded according to a cut-off of 3 standard deviations (i.e., z = |3|).Note that we also changed the threshold from 2.5 to 3.5 to confirm that the results were not dependent upon the choice of the threshold (Supplementary materials 1-3).To be sure that we were actually capturing the critical dynamics (if present), we binned the time series.Similarly to Sorrentino et al. (2021a), we estimated the suitable time bin length by computing the branching ratio σ for each subject, each neuronal avalanches, and each time bin duration (Haldeman and Beggs, 2005).Indeed, a branching ratio 1 typically indicates a critical process.Specifically, the branching ratio was calculated as: where σ 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, and n events (j) is the total number of events in the j-th bin.We then (geometrically) averaged the results over all avalanches as follows: However, we repeated our analysis exploring different time bins, ranging from 1 to 5 obtaining similar results (see Supplementary materials 2-3).For each avalanche, an avalanche pattern was defined as the set of all brain regions that were above threshold.

Functional repertoire
The functional repertoire represents the number of unique avalanche patterns expressed during the recording (Sorrentino et al., 2021a).Unique means that each avalanche pattern is counted only once to define the avalanche size of the functional repertoire (i.e., repetitions are discarded).

Transition matrices
We built an ATM in which its element (i, j) represented the probability that region j was active at time t + δ, given that region i was active at time t.That is, 1 ATM was build per each avalanche, and then they were averaged element-wise to obtain a subject-specific ATM.The ATMs were then averaged per participant and symmetrized.Please note that ATMs are extremely sparse matrices.In the subject-specific ATM, each ijth entry is computed from those avalanche-specific ATMs when the ij-th edge was recruited, and not from the others.As such, the quantity of data that contributes to the estimate is edge-specific.However, since avalanches are rare, fast-lived events, each edge typically refers to 5%-10% of the data points, sampled in a in homogeneous fashion (nonperiodically) since the sampling is done according to when the avalanches were occurring (for more details, see Sorrentino et al., 2021c).We explored the frequency-specific ATMs, by filtering the source-reconstructed signal in the classical frequency bands (delta, 0.5-4 Hz; theta, 4-8 Hz; alpha, 8-13 Hz; beta, 13-30 Hz; gamma, 30-48 Hz) since recent lines of evidence suggest a coexistence between the frequencyspecific oscillatory activity and large-scale aperiodic perturbations.These 2 kinds of activities (i.e., periodic and aperiodic) operate simultaneously, are embedded into each other, and affect each other (Beggs and Plenz, 2003;di Santo et al., 2018;Gireesh and Plenz, 2008;Lombardi et al., 2023;Poil et al., 2012).However, we wish to stress that we filtered the data in the classical frequency bands since their physiological significance is well established.However, we cannot rule out that the filtering introduces distortions on the aperiodic activities.
Thereafter, based on the frequency-specific ATMs, we estimated the topology of the avalanches by calculating the minimum spanning tree (MST), a spanning tree (a subgraph that, while connecting every node, has no cycles) whose edges have the least total weight.In particular, we calculated the MST using Kruskal's algorithm (Kruskal, 1956) and extracted both global and nodal topological parameters.The former were the degree divergence (DD), which represents the amplitude of the degree distribution (i.e., how much the most central nodes differ, in terms of number of connections, from the least central nodes); the leaf fraction (LF), which is defined as the fraction of nodes with degree = 1; the tree hierarchy, which conveys the trade-off between an efficiently connected network and its resiliency to targeted attacks, and the assortativity, which represent the tendency of a node to be connected to similar nodes (in terms of degree) in a given network (Jacini et al., 2018;Stam et al., 2014;Tewarie et al., 2015).Then, as nodal parameters, we calculated the betweenness centrality and the degree, which represent the topological importance of a node within the network (Rubinov and Sporns, 2010) ( Fig. 1).

Multilinear regression analysis
We hypothesized that the alterations of brain dynamics could be used to predict the patients' clinical impairment.Therefore, we built a multilinear regression model in which the neuropsychological tests represented the dependent variable while age, gender, education, and brain dynamics features were the predictors.Multicollinearity was assessed through the variance inflation factor.To validate our approach, we performed k-fold cross-validation, with k = 5 (Varoquaux et al., 2017).Specifically, k iterations were performed to train our model, and at each iteration, the k-th subgroup was used as a test set.

Statistical analysis
Statistical analysis was carried out in MATLAB 2021a.The comparison between patients and HC was performed through permutation testing.Data labels of both groups were shuffled 10,000 times.At each iteration, 2 random groups were generated, and the absolute difference between the means of the 2 groups was calculated.Hence, we obtained a null distribution of 10,000 random differences, to which we compared the observed absolute difference between the means of patients and HC.Finally, a possible relationship between the topological features of brain dynamics and the clinical conditions (hippocampal volumes and cognitive impairment) of the patients was explored through Spearman's correlation.The results were corrected for multiple comparisons across both frequency bands and parameters using the false discovery rate (FDR) method (Benjamini and Hochberg, 1995) and the significance level was set at p value < 0.05 (corrected).Finally, we have used phase-surrogates in order to demonstrate that the changes in topology genuinely depend upon the interactions between regions and are not to be expected by the regional dynamics alone.Similarly to Sorrentino et al. (2023b), we generated surrogates that disrupt the dependencies between regions, while preserving the local (i.e., regional) activities.In short, starting from the real, broadband, source-reconstructed time series, we have performed a Fourier transform and shifted the phase of each frequency component by a random amount (while preserving the Hessian symmetry of the power spectra), and then antitransformed the signal.This way, we have obtained surrogate real signals which are characterized by the same power spectra as the original data, while the relationships between frequencies are disrupted.After doing this, we performed the same analyses that we had previously reported on the surrogate data.

Cohort characteristics
As shown in Table 1, the comparison between the demographic characteristics of the MCI and the HC groups revealed no significant differences in both age and education level.With respect to the structural information, only 30 MCI patients out of 32 and 22 out of 32 HC completed the MRI scan.MCI patients exhibited a reduction of hippocampal volumes compared to the HC (p < 0.001).As expected, MCI displayed a worse cognitive performance.Specifically, the comparison between the neuropsychological evaluation of the 2 groups showed a significant difference in the: Mini Mental State Examination (p = 0.012), Rey's 15-word immediate recall (p < 0.01), Rey's 15-word delayed recall (p < 0.01), word fluency (p < 0.01), Raven's 47 progressive matrices (p < 0.01), immediate visual memory (p < 0.01), and the FCSRT (immediate free recall, immediate total recall, delayed free recall, delayed total recall, and index of sensitivity of cueing) (p < 0.001).No significant differences were found in the evaluation of the Beck Depression Inventory, Frontal Assessment Battery, Phrase construction, Freehand copying of drawings, and Constructive apraxia with landmarks (see Table 2).

Functional repertoire
In order to analyze the size of the functional repertoire (i.e., the number of unique and different avalanches configurations in each participant), we used source-reconstructed MEG data.As shown in Fig. 2, the MCI group showed a reduction in the number of unique patterns in the theta band (pFDR = 0.047), revealing a restricted functional repertoire as compared to the HC.These analyses have been repeated on the surrogate data (Supplementary materials 4), and the differences disappear.

Transition matrices topology
To study avalanche topology, we built an ATM (Sorrentino et al., 2021c) for each subject.Then, we computed the MST for each ATMs (Wijk et al., 2010).Patients showed higher values of LF and DD compared to the healthy participants (Fig. 3).Specifically, with respect to the LF, we found a statistically significant difference between the 2 groups in the delta (pFDR = 0.020) and theta bands (pFDR = 0.041).Similarly, we found a statistically significant difference of the k between the 2 groups in the delta (pFDR = 0.006), theta (pFDR = 0.046), and gamma bands (pFDR = 0.020).We did not find any statistical difference neither for the nodal parameters (i.e., betweenness centrality and degree) nor for the remaining global parameter (i.e., the tree hierarchy).Again, these results are not present in the surrogate data (Supplementary material 4).

Correlations with clinical parameters
Spearman's correlation was used to determine whether there was a relationship between topological network features and both the episodic memory (assessed through the FCSRT) and the structural parameters (i.e., hippocampal volumes).We focused on the FCSRT due to its sensitivity in detecting the hippocampal damage and amnesic deficits even in the earliest phase of disease (Auriacombe et al., 2010;Dubois et al., 2014).Firstly, we found a significant positive correlation, after FDR correction, between the DD in the theta band and the FCSRT immediate total recall (r = 0.536; pFDR = 0.046).In addition, the DD in the theta band positively correlated with the right hippocampal volume (r = 0.438; pFDR = 0.030) (Fig. 4).Lastly, we found a significant positive correlation between the FCSRT immediate total recall and both right (r = 0.553; pFDR = 0.001) and left (r = 0.463; pFDR = 0.009) hippocampal volumes (data not shown).No significant correlation was found between the topological parameters and the other neuropsychological test scores and the left hippocampal volume.

Predictive model: multilinear regression analysis
Taken into consideration the relationship between the DD (i.e., the width of the degree distribution) and both the memory impairment and the right hippocampal atrophy, we asked ourselves whether the width of the degree distribution could improve the prediction of the clinical impairment (as measured via the FCSRT immediate total recall).To this end, we built a multilinear regression model and we validated it using a k-fold cross-validation approach.The model also contained 4 nuisance predictors, that is, age, education level, gender, and the right hippocampal volume (see Fig. 5).We added the right hippocampal volume as a predictor given its role in memory and its involvement in MCI and AD (Apostolova et al., 2006;Devanand et al., 2012;Sarazin et al., 2010;Yavuz et al., 2007;Zammit et al., 2017).Our aim was to observe whether the DD improved the prediction over the FCSRT immediate total recall as compared to using demographics and structural changes (i.e., right hippocampal volume) alone.We found that the DD in the theta band significantly predicted the FCSRT immediate total recall (p = 0.041; R 2 = 0.47; β = 2.7).Furthermore, the only other significant predictor was the education level (p = 0.012; β = 0.448) (Fig. 5).As an extra check, we rebuilt the same predictive model, but this time we switched the order of the last 2 predictors (i.e., the DD and the right hippocampal volume).Our aim was to visually assess the residual predictive power of the hippocampus after the DD was added.As shown in the supplementary materials (Fig. S5), the hippocampal volume did not significantly contribute to the prediction of the functional impairment after the functional DD was taken into account.

Discussion
In the present work, we used source-reconstructed MEG data in a cohort of 32 MCI patients and 32 healthy participants to explore the topological organization of the large-scale brain dynamics and its relationship with the clinical impairment.Firstly, we hypothesized that MCI participants would show more stereotyped brain dynamics and, consequently, such stereotyped brain dynamics might result in a reorganization of the spreading of the large-scale perturbations, as captured by changes in the functional topology.We used neuronal avalanches to quantify the flexibility of fast brain dynamics (Sorrentino et al., 2021a).The brain flexibility is given by its ability to generate a large number of functional configurations, that likely conveys a large number of different ways by which brain regions can interact among themselves (Shine et al., 2016;Zalesky et al., 2014).These large-scale interactions manifest themselves as bursty activity (i.e., neuronal avalanches), and the number of different ways neuronal avalanches propagate across the brain defines the size of the functional repertoire (Chialvo, 2010;Sorrentino et al., 2021a).
Our results showed that MCI patients exhibited a restricted functional repertoire (i.e., lower number of unique avalanche patterns) in the theta band compared to HC.Thus, MCI patients displayed a reduction of their brain flexibility, which implies less reconfigurations of brain activity and a more stereotyped brain dynamics (Sorrentino et al., 2021a).This is in line with previous reports that used the same approach in other neurodegenerative diseases such as Parkinson's disease (Sorrentino et al., 2021a) and amyotrophic lateral sclerosis (Polverino et al., 2022), suggesting that the loss of flexibility might be a common dynamical change induced by different neurodegenerative processes.
Of note, our results were evident mainly in the theta band and were undetectable from the broadband signals (Supplementary materials 6).We speculate that the aperiodic bursts, which are the effect of the resonances occurring among cross-regional frequencyspecific activities, possess themselves an emergent characteristic  time-scale, which makes them detectable after appropriate filtering.However, further studies will have to systematically investigate the rhythmicity of the aperiodic activities and the effect of different filtering strategies.
We performed surrogate data analysis to partially test this hypothesis by disrupting the cross-regional resonances while keeping the regional power spectra unchanged.In other words, we disrupted the (phase) relationships among different frequency bins.We performed this procedure over the entire frequency spectrum, without selecting any specific frequency a priori.When we repeated the whole analysis pipeline on surrogate data, including the filtering procedure, we show that the differences between MCI patients and controls are no longer evident in the case of surrogates.We conclude that the resonances generated by all the frequency-specific activities produce higher-order bursts that spread with characteristic timescales.On a speculative note, this might be due to specific subpopulations that operate at particular timescales, which is reminiscent of the well-known bursty dynamics occurring in Parkinson's disease in the beta-frequency range (Bonaiuto et al., 2021).
Then, we generated the ATMs which provide information about the spatio-temporal avalanches dynamics estimating the probability that 2 brain regions will move away from their baseline activity in 2 consecutive time frames (Sorrentino et al., 2021c).By applying the MST on the ATMs, we proceeded to investigate their topological properties.The comparison between MCI patients and HC revealed that the former exhibited higher values of both LF (in the delta and theta bands) and DD (in the theta and gamma bands) with respect to the latter.Higher LF values (i.e., fraction of nodes with degree = 1) are indicative of a more integrated network (as in a star-like topology), suggesting the idea that, when the LF increases in the patients, the network shifts toward a more centralized organization (Boersma et al., 2013).The DD (i.e., the width of degree distribution) is a measure of the network's synchronisability and its resilience against pathological events (Tewarie et al., 2015).Higher values of DD convey the presence of high-degree nodes (Boersma et al., 2013), which, again, is consistent with a more centralized topology of the network.In line with previous evidence (Sorrentino et al., 2018;Stam, 2014), one might speculate that the higher centrality might be due to hubs compensating for the impairment of more peripheral nodes.As a consequence, the communication between brain regions becomes less efficient.Hence, MCI participants can put in place more effective compensation mechanisms, whereby the reduction of lower-degree nodes might be effectively compensated by the more central hubs (Rucco et al., 2019).A possible unifying explanation of the current results might lay in a functional compensatory mechanism which would be in place before the structural changes.In other words, one could speculate that, in the earliest phase of a neurodegenerative disease, functional anomalies, for example, as in the case of Aβ-mediated synaptic dysfunctions in the initial stages of AD (Zhang et al., 2022), provoke a whole-brain functional network rearrangement in order to maintain proper cognitive performances.Conversely, as the disease progresses, this "functional" failure is flanked by a progressive loss of gray matter that becomes preponderant in the advanced stage of the disease, when multiple brain areas are involved (Cipriano et al., 2022).What stated so far is in line with our results showing that an increase in DD reflects a better cognitive performance and is also correlated with greater right hippocampal volume.Indeed, the higher the DD, the better the FCSRT immediate total recall scores and the larger the hippocampal volume.In addition, we found a positive correlation between the FCSRT immediate total recall scores and the right and left hippocampal volumes.Our results are in line with previous evidence (Jacini et al., 2018;Wang et al., 2006), which demonstrated the relationship between the loss of memory efficiency and hippocampal atrophy.For instance, a longitudinal study conducted by Stoub et al. (2010) showed that, over 5 years, MCI patients showed a progressive cognitive impairment directly related to changes of hippocampal volume.
Specifically, our results revealed a correlation between the topology (i.e., DD in theta band) of the MCI and the right hippocampal volume.According to Cai and colleagues, MCI patients show a reduction of their FC in the right hippocampus (which is also characterized by a greater rate of atrophy with respect to the left one).Importantly, these alterations may be related to a worse cognitive performance expressed by a deficit of the episodic memory retrieval Fig. 5. Clinical impairment prediction in the theta band.Multilinear regression analysis with k-fold cross validation was performed to verify the ability of the degree divergence to predict the memory impairment assessed by the FCSRT immediate total recall.The left column displays the explained variance obtained by adding the predictors (age, education level, gender, right hippocampal volume, and the degree divergence in the theta band).The significant predictors are highlighted in bold while the significant p-value is indicated with *(p < 0.05), **(p < 0.01).The central column displays the comparison between the predicted and the actual values of the responsive variable validated through the k-fold cross validation.Lastly, the right column shows the distribution of residuals, which represent the standardization of the difference between the actual and predicted values for the FCSRT immediate total recall.Abbreviations: FCSRT, free and cued selective reminding test.(Cai et al., 2017).This is in agreement with Heckers et al. who showed that the best retrieval of encoded words through a semantic cue was associated with higher right hippocampal activation (Heckers et al., 2002).However, we would like to specify that the brain network alterations in MCI are not constrained to the right hippocampus, but extend to the left hemisphere, reflecting, similarly to the right hemisphere, poorer cognitive performance (Bai et al., 2009).In this regard, O'Neill et al., in a recent study, analyzed the source reconstruction in MEG with regard to hippocampal activity, finding a possible bias of functioning with regard to activity reconstruction performed with LCMV beamformer.In particular, it is observed that the LCMV beamformer, in cases of correlated source suppression (e.g., bilateral activation of hippocampus), may not be able to reconstruct the entire activity (O'Neill et al., 2021).These limitations warrant caution in interpreting the results.In particular, we cannot rule out that these technical limitations might have a role in the fact that we did not find any functional differences in the analysis of hippocampal dynamics, and that we only found a unilateral correlation between DD and hippocampal volume.
In the current paper, we wished to take advantage of the information about the fast brain dynamics, in order to further improve the individual predictions of clinical impairment (Bosboom et al., 2006;Chen et al., 2021).To this end, we built a multilinear regression model which demonstrated the predictive power of the DD in the theta band over the FCSRT immediate total recall.Including the DD in the model improved predictions over other nuisance predictors such as age, gender, education level, and the right hippocampal volume.Interestingly, as shown in supplementary materials (Fig. S5), the structural information alone is unable to predict the clinical impairment (as assessed by the FCSRT immediate total recall), while functional rearrangement (defined by the higher DD) carries predictive power.
In summary, these results show that the functional data provide a more comprehensive information, which also takes into account potential compensatory mechanisms which might be overlooked when focusing purely on structural data (i.e., atrophy) (Cipriano et al., 2023).These findings show that dynamical features predict functional abilities (i.e., behavior).Besides the DD, only the education level was a significant predictor for the FCSRT immediate total recall.Although the FCSRT scores are corrected by taking into account the education level, it is quite obvious that such a relationship would exist since it is expected that highly educated participants are also those who exhibit better cognitive performance (Frasson et al., 2011).Finally, it is worth mentioning that the topological properties of both MCI patients and healthy participants were obtained from the ATMs, that characterize the spatio-temporal evolution of neuronal avalanches.In other words, our evidence corroborates the clinical and behavioral relevance of aperiodic activity at the fast timescales.

Conclusion
In the current work, we applied the novel framework based on ATMs to investigate the spatio-temporal dynamics of neuronal avalanches in MCI and HC.As we hypothesized, MCI patients displayed a reduction of their brain dynamic flexibility, which resulted in suboptimal functional topology and predicted the memory impairment.Based on these findings, we hope that the proposed framework may be helpful in monitoring the development of the disease and effectively capturing subtle information derived from dynamical analysis to achieve a more thorough phenotypical assessment of patients.

Fig. 1 .
Fig. 1.Pipeline overview.(A) 1. Registration of neuronal activity through MEG; 2. Raw sensor signals including physiological artifacts; 3. Cleaned sensor signals (without physiological artifacts); 4. MRI of the subject; 5. Beamforming obtained from the coregistration of MRI and MEG signals.(B) Source-reconstructed temporal series.The blue lines represent the z-score activity of a brain region, and the light blue embossed boxes represent the time frame in which a neuronal avalanche occurred.Specifically, the red dots indicate the frame in which the time series was above threshold (z-score > 3).(C) Schematic representation of a neuronal avalanche.A neuronal avalanche occurred when at least 1 region was above threshold and ended when all the brain regions returned below threshold.The avalanche time-evolution matrix displays which regions (only 4 are displayed for visualization purposes) were above threshold (blue boxes) during each specific time frame.Collapsing the time outlines an avalanche pattern that includes regions that were active for at least 1 time frame during the avalanche.As a whole, the patterns represent the functional repertoire, a measure of brain flexibility.(D) Schematic representation of a transition matrix.The light blue squares indicate that the region i was above threshold 3 times during the avalanche.Region j was active, after the activation of region i, only in 2 cases out of the 3 taken into consideration (as indicated by the green arrows).Thus, the probability that region j was active after the activation of region i was 2/3.Averaging across avalanche-specific transition matrices creates the avalanche transition matrix of an individual.Finally, a topological analysis was conducted to highlight the network properties of the brain activity propagation.Abbreviations: MEG, magnetoencephalography; MRI, magnetic resonance image; MST, minimum spanning tree.

Fig. 2 .
Fig. 2. Functional repertoire comparison.The violin plots represent the comparison of the number of unique avalanche patterns between HC and MCI patients in theta band.MCI patients display a reduced functional repertoire compared to the HC.The median is represented by the white dot, the gray bar in the middle of the violins shows the first and the third quartile, and the thin gray line is representative of the 95% confidence interval.The outliers are individually represented by the single dots.*p < 0.05.Abbreviations: HC, healthy controls; MCI, mild cognitive impairment.

Fig. 3 .
Fig. 3. Topological parameters comparison.Violin plots represent the comparison between HC and MCI patients for topological parameters obtained from the transition matrices.MCI patients display higher values of leaf fraction in delta and theta bands with respect to HC as well as higher values of degree divergence in delta, theta, and gamma bands.*p < 0.05; **p < 0.01.Abbreviations: HC, healthy controls; MCI, mild cognitive impairment.

Fig. 4 .
Fig. 4. Correlation between MCI clinical features and degree divergence.(A) Spearman's correlation between the DD and the FCSRT immediate total recall.The positive correlation coefficient shows that, as the DD increases, the FCSRT score increases as well.(B) Positive correlation between the DD and the right hippocampal volume.Higher values of DD correspond to a greater preservation of the right hippocampal volume of the patients.For panel B, the Spearman's correlation analysis was performed taking into account the only 30 MCI patients who complete the MRI scan.(C) For completeness, we reported the not significant correlation between the left hippocampal volume and the DD in theta band.Abbreviations: FCSRT, free and cued selective reminding test; FDR, false discovery rate.
The structural information was available for 22 out of 32 healthy controls and for 30 out of 32 MCI patients.