Brain functional networks become more connected as amyotrophic lateral sclerosis progresses: a source level magnetoencephalographic study.

This study hypothesizes that the brain shows hyper connectedness as amyotrophic lateral sclerosis (ALS) progresses. 54 patients (classified as "early stage" or "advanced stage") and 25 controls underwent magnetoencephalography and MRI recordings. The activity of the brain areas was reconstructed, and the synchronization between them was estimated in the classical frequency bands using the phase lag index. Brain topological metrics such as the leaf fraction (number of nodes with degree of 1), the degree divergence (a measure of the scale-freeness) and the degree correlation (a measure of disassortativity) were estimated. Betweenness centrality was used to estimate the centrality of the brain areas. In all frequency bands, it was evident that, the more advanced the disease, the more connected, scale-free and disassortative the brain networks. No differences were evident in specific brain areas. Such modified brain topology is sub-optimal as compared to controls. Within this framework, our study shows that brain networks become more connected according to disease staging in ALS patients.

Introduction

Recently, autoptic, neuropsychological, genetic and neuroimaging evidence has suggested that amyotrophic lateral sclerosis (ALS) involves regions beyond the primary motor cortex (Turner and Swash, 2015).

According to this evidence, fMRI has identified wide areas of both increased and decreased connectivity across the whole brain (Agosta et al., 2011, Agosta et al., 2013; Chiò et al., 2014). Several hypotheses have been taken into account to explain such evidence, including the interpretation of hyper-connectivity as a compensation mechanism or as a by-product of the loss of inhibitory circuitry (Turner and Kiernan, 2012). Furthermore, it was proposed that damage to peripheral areas of the brain might lead more central areas to process information that can no longer be handled locally, resulting in the overload of such areas (Stam, 2014). In primary lateral sclerosis, it was reported that patients with the greatest clinical disability also displayed the highest functional connectivity (Agosta et al., 2014). Interestingly, such increased functional connectivity was widespread, since it was present in the sensorimotor, premotor, prefrontal and talamic regions (Agosta et al., 2014). Furthermore, it was shown in ALS, with combined structural and functional MRI, that regions with lower structural connectivity displayed increased functional connectivity (Douaud et al., 2011). Moreover, patients with a slow rate of progression displayed lower connectivity, thus being more similar to the controls (Douaud et al., 2011). Recently, in a small, MRI based study, cortical excitability (evaluated by the threshold-tracking transcranial magnetic stimulation paradigm) negatively correlated with functional connectivity (Geevasinga et al., 2017). Furthermore, Verstraete et al., by exploring abnormalities of MRI structural and functional connectivity in a cohort of ALS patients compared to healthy controls, noted that the connectedness of the network related directly with the progression rate (Verstraete et al., 2010). In a large, fMRI based study, Schulthess et al. revealed increased functional connectivity, with a topography linked to the spreading of the pTDP-43 pathology (Schulthess et al., 2016). The reported evidence might be compatible with the idea that hyper connectedness might relate to neuronal damage in ALS (Trojsi et al., 2017).

Beside fMRI, relevant information on the functioning of brain networks can be obtained using neurophysiological techniques, such as electroencephalography (EEG) and magnetoencephalography (MEG), since they directly capture the electrical/magnetic activity of the neuronal ensembles (Lopes Da Silva, 2013). In ALS, increased EEG coherence in the theta and gamma band was shown to correlate with structural MRI changes (Nasseroleslami et al., 2017) and EEG based network topological metrics related to disease burden (Fraschini et al., 2016b).

MEG systems measure magnetic fields produced by neuronal activity. Such fields are minimally distorted by the layers surrounding the brain, allowing for a temporally and spatially accurate reconstruction of the neural signals within the brain (source space) (Baillet, 2017).

Recently, a source level MEG study based on power confirmed that hyperactivation might be a relevant feature in ALS during eyes opened resting state (Proudfoot et al., 2018). Such evidence would confirm further that the hyperactivation would relate to neuronal damage.

However, MEG can be further exploited extracting the phases of the signals in order to evaluate the amount of information exchanged between brain areas. More precisely, the estimation of the phases allows the quantification of true synchrony, defined as a constant phase difference between time-series in isofrequency (Tass et al., 1998). We chose the phase lag index (PLI) (Stam et al., 2007) to estimate functional connectivity between brain areas, since the PLI quantifies synchronization between time series while being entirely unaffected by the amplitude of the signal. The PLI was also chosen since it is insensitive to volume conduction (at the cost of discarding true zero-lag interactions). Interestingly, the fact that the PLI estimates synchronization rather than simultaneous activations of brain areas, makes the information provided in this paper complementary to the fMRI based evidence and power-based MEG analyses.

Some of the properties of the interactions among brain areas can be captured using Graph Theory (Bullmore and Sporns, 2009). In fact, the human brain can be modelled as a network, with the brain areas as nodes and their interactions as edges. However, comparing such metrics is non-trivial, as they are influenced by network size, thresholding or edge density, not allowing for a purely topological interpretation of the results (van Wijk et al., 2010). The minimum spanning tree (MST) algorithm allows the computation of statistically comparable metrics, while retaining most of the information about the original network (Tewarie et al., 2015). The combination of the PLI and the MST, while it is on the one hand costly in the sense that it might discard some information, on the other hand reduces to a minimum the detection of either artefactual or uninterpretable differences between groups.

On the basis of the available evidence, we hypothesized that the topological alterations associated to different stages of ALS may involve vast areas of the brain and not be confined to motor areas. Secondly, given the evidence showing that hyper-connectedness is related to atrophy and disease progression, we hypothesize that, as the patients reach more advanced stages, the brain network will show a more connected topology accordingly. A more integrated topology implies that it is quicker and/or less costly, on average, to move among the nodes of the network. Lastly, we hypothesize that the alterations will be spread across frequency bands. To test our hypothesis, a large cohort of ALS patients and healthy controls underwent clinical evaluation, MRI and MEG scan. The population was classified into “early” and “advanced” stage on the base of the King's disease staging system (Balendra et al., 2015). In particular, considering the current unavailability of validated markers of disease progression, we chose to base our analyses on such clinical staging system since it seems to reproduce the curvilinear course of disease progression typical of ALS (Gordon et al., 2010) especially in the later stages of the disease (Balendra et al., 2015). Based on the PLI and the MST, we computed topological metrics, focusing on features such as scale-freeness, assortativity and connectedness of the brain networks, as well as on the betweenness centrality of each area, in order to test if such parameters would differ among controls and ALS patients at different disease stages.

Material and methods

Participants

Fifty-four right-handed and native Italian speakers patients (39 males, 15 females; mean age ± SD, 58.84 ± 12.14) with probable or definite ALS, according to the revised El-Escorial criteria of ALS (Brooks et al., 2000), were consecutively recruited at the ALS Center of the First Division of Neurology of the University of Campania “Luigi Vanvitelli” (Naples, Italy). The patients were classified according to the King's disease staging system (Balendra et al., 2015) that is based on the appearance of sequential clinical milestones during the ALS course and does not include cognitive information (i.e., stage 1 = impairment of one body site; stage 2 = impairment of two body sites; stage 3 = impairment of three body sites; stage 4 = non-invasive ventilation or percutaneous endoscopic gastrostomy). In order to improve the group sizes, we classified as “early stage” the patients belonging to King's stage 1 and 2 (9 and 17 patients, respectively) and as “advanced stage” the patients in King's stages 3 and 4 (12 patients in each group). None of the patients showed any mutation in any of the following genes: SOD1, TARDBP, FUS/TLS and C9ORF72. More clinical details of participants included in the final analysis (4 patients were dropped out of the analysis as they did not have enough high quality MEG data) are given in Table 1. Twenty-five controls (16 males and 9 females; mean age ± SD, 57.00 ± 9.35), enrolled by “word of mouth”, were age-, −gender and education-matched with the ALS patients. To be included in this study, all participants had to satisfy the following criteria: a) to have no major medical illnesses and not to abuse substances or use medication that could interfere with MEG/EEG signals; b) to show no other major systemic, psychiatric, or neurological illnesses; and c) to have no causes of focal or diffuse brain damage, including lacunae, and extensive cerebrovascular disorders at routine MRI. Four patients did not have a MEG recording of sufficient quality and were discarded from the study. The study protocol was approved by the local Ethics Committee. All participants gave written informed consent.

Table 1

Detailed characteristic of patients and controls used for the analysis.

Parameters	ALS “advanced” patients mean (SD) (n = 24)	ALS “early” patients mean (SD) (n = 26)	Controls mean (SD) (n = 25)
Demographic and clinical measures
Age	59.96 (13.89)	57.50 (10.76)	57 (9.35)
Male/Female	19/5	18/8	16/9
Education	10.46 (4.51)	10.19 (4.09)	11 (4)
Disease duration (months)	61.83 (60.23)	28.77 (20.69)
ALSFRS-R score	30.70 (8.79)	41.15 (4.80)
UMN score	8.22 (5.13)	6.46 (4.58)
Site of onset	6 bulbar	5 bulbar
	6 UL	11 UL
	9 LL	9 LL
	2 UL and LL	1 UL and LL
	1 respiratory	0 respiratory
Phenotype	9 classic	9 classic
	6 predominant LMN	12 predominant LMN
	9 predominant UMN	5 predominant UMN
Neuropsychological parameters
ECAS test (total score)	83.13 (28.08)	91.50 (21.09)

ALSFRS-R = Amyotrophic Lateral Sclerosis Functional Rating Scale-Revised; ECAS = Edinburgh Cognitive and Behavioural ALS Screen; LL = Lower Limb; LMN = Lower Motor Neuron; UL = Upper Limb; UMN = Upper Motor Neuron.

MRI acquisition

For 41 patients and 22 controls, 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 the MEG recording or a minimum of 21 days earlier (within one month). In particular, 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. The remaining participants did not complete the MRI because of the difficulty to lie down or refused to perform the MRI scan. In this case, we used a standard template MRI.

MEG acquisition

Subjects underwent a magnetoencephalographic examination in a 163-magnetometers MEG system (Rombetto et al., 2014) placed in a magnetically shielded room (AtB Biomag UG - Ulm - Germany). The position of four position coils and of four reference points (nasion, right and left pre-auricolar and apex) were digitized before acquisition. The brain activity was recorded for 10 min, alternating between eyes closed and eyes opened condition every 2.5 min, so as to minimize the chances of drowsiness. The instruction to close/open eyes was given using a standardized sentence communicated through intercom. The head position was recorded at the start of each segment of the recording. The order between eyes closed and eyes opened segments was randomized for both patients and controls. After the anti-aliasing filter, the data were sampled at 1024 Hz. Subsequently, a 4th order Butterworth IIR band-pass filter has been applied in order to remove components below 0.5 and beyond 48 Hz. The filter is implemented within the Fieldtrip toolbox (Oostenveld et al., 2011). During the acquisitions, electrocardiogram (ECG) and electrooculogram (EOG) were also recorded (Gross et al., 2013).

Preprocessing

Eyes-closed trials were selected for further analysis. Principal component analysis (PCA) was performed to reduce the environmental noise (Sadasivan, 1996). Specifically, the filter orthogonalizes the reference signals to obtain a basis and projecting the brain sensors on the basis of the noise, removing the projections to clean the data (de Cheveigné and Simon, 2007). We adopted the PCA filtering implementation available within the Fieldtrip Toolbox (Oostenveld et al., 2011). Subsequently, noisy channels were removed manually through visual inspection of the whole dataset by an experienced rater (Gross et al., 2013). On average, 130 ± 2 channels have been used. For each subject, supervised independent component analysis (ICA) (Barbati et al., 2004) was performed to eliminate the ECG (generally one, rarely two components) and the EOG (none, rarely one components) component from the MEG signals. It is important to notice that ICA modifies to some extent the phases of the signal (Bridwell et al., 2018; Montefusco-Siegmund et al., 2013; Shou and Ding, 2014). In order to test if this would affect our results, one population (i.e. controls) was also cleaned by visual inspection without ICA. The results have been then compared to those obtained after ICA had been applied, and no difference was evident either at visual inspection or after statistical analysis (see Supplementary material 1 for details).Ten epochs of 8 s for each subject that did not contain artefacts (either system related or physiological) or strong environmental noise were selected. As stated earlier, 4 patients did not meet the criteria, and 50 patients went on for further analysis. The length of 8 s is a trade-off between the need to have enough cleaned epochs, to avoid drowsiness (Gross et al., 2013) and to obtain a reliable estimate of the connectivity measure (Fraschini et al., 2016a). However, in order to make our results more robust, also 6 s long and 10 s long epochs have been analyzed and no difference was observed (data not shown).

Source reconstruction

All the processing related to the beamforming procedure has been done using the Fieldtrip toolbox (Oostenveld et al., 2011). Initially, the subject's fiducial points were used to coregister the MEG acquisition to the native MRI of the subjects by computing the transformation matrix. In case the native MRI was not available, a template was adopted. Subsequently, the volume conduction model proposed by Nolte (Nolte, 2003) was considered and the Linearly Constrained Minimum Variance (LCMV) beamformer (Van Veen et al., 1997) was implemented to reconstruct time series related to the centroids of 116 regions-of-interest (ROIs), derived from the Automated Anatomical Labeling (AAL) atlas (Gong et al., 2009; Hillebrand et al., 2016; Tzourio-Mazoyer et al., 2002). We considered only the first 90 ROIs, excluding those corresponding to the cerebellum given the low reliability of the reconstructed signal in those areas. For each source, we projected the time series along the dipole direction that explains most variance by means of singular value decomposition (SVD). Source space time series were re-sampled at 512 Hz. For each subject we selected by visual inspection at least ten epochs of 8 s free from artefacts (both system related or physiological) or excessive environmental noise. An overview of the whole pipeline is provided in Fig. 1.

Fig. 1

Data analysis pipeline.A.Neuronal activity in the sensor space as recorded by magnetoencephalography (MEG). Alpha activity has been represented in the image as an example. B. Magnetic resonance (MR) of the subject. C. The MR and the MEG sensors are coregistered (i.e. they are in the same space). D. A model of the brain volume is created. E. The Activity on the sensor level is back projected onto brain regions (based on the AAL atlas). F. The phase lag index (PLI) is estimated in a pairwise fashion between each pair of 90 (AAL-based) brain regions. G. Based on the Kruskal's algorithm, the minimum spanning tree is reconstructed. H. The minimum spanning tree is represented graphically, where raws and columns of the matrix in G are represented as red dots, and the entries of the matrix in G are represented as lines. Once a frequency specific MST has been obtained, it will be used to compute topological metrics.

Connectivity analysis

The phase lag index (PLI) was chosen to measure synchrony between brain areas (Stam et al., 2007). The PLI was computed using BrainWave software [CJS, version 09.152.1.23, available from http://home.kpn.nl/stam7883/brainwave.html]. The instantaneous phases were computed using a Hilbert transform of the time series, and the PLI was used to estimate the asymmetry of the distribution of the phase differences (Δφt) between the two time series as:where ‘<.>’ indicates the mean value, ‘sign’ stands for the signum function, ‘|.|’ denotes the absolute value and ‘t’ are the discrete time-steps. A PLI of 0 indicates completely symmetric distribution of the phase differences, or phase differences of zero (mod π), and 1 indicates perfectly asymmetric distribution of the phase differences (Stam et al., 2007). This metric has been shown to be insensitive to volume conduction (Stam et al., 2007), since zero lag interactions do not contribute to the estimate of the PLI. However, this implies that true zero lag interaction will also be neglected. By computing PLI for each couple of brain regions, we obtained a 90 × 90 weighted adjacency matrix for each epoch and for each subject, in all of the frequency bands: delta (0.5–4 Hz), theta (4.0–8.0 Hz), alpha (8.0–13.0 Hz), beta (13.0–30.0 Hz) and gamma (30.0–48.0 Hz).

Network analysis

The network analyses were performed using BrainWave software as before. The weighted adjacency matrix was used to reconstruct a network or complete weighted graph, where the 90 areas of the AAL atlas are represented as nodes, and the PLI values form the weighted edges. A frequency specific minimum spanning tree was calculated for each epoch. Since we were interested in the strongest connections, for the construction of the MST the edge weight was defined as 1/PLI. In fact, Kruskal's algorithm first ranks the links in ascending order and then constructs the network by adding one link at a time, discarding links that would form a loop. The algorithm proceeds until all nodes are connected resulting in a loop-less graph with N nodes and M = N-1 links.

We used the minimum spanning tree (MST) to avoid some of the biases in traditional network analyses (Stam, 2014; Stam et al., 2014; Tewarie et al., 2015; van Wijk et al., 2010), to allow for an unbiased topological interpretation of the results (Stam et al., 2014). Based on the MST, we calculated the leaf fraction, the degree divergence and the degree correlation and the tree hierarchy. The leaf fraction is defined as the fraction of nodes with degree of 1 (Boersma et al., 2013), providing an indication of the integration of the network. A higher leaf fraction implies that the network tends toward a star-like topology, where the nodes are on average closer to each other as compared to a more line-like topology (see Fig. 3). The degree divergence is a measure of the broadness of the degree distribution, related to the resilience against attacks (Tewarie et al., 2015). The degree correlation measures to what extent the degree of a node is influenced by the degree of its neighbors. Degree correlation is computed using the Pearson correlation coefficient of the degrees of pairs of connected nodes. If the correlation is positive the graph is called assortative, and if the correlation is negative the graph is called disassortative (Boersma et al., 2013; Newman, 2002). Finally, the tree hierarchy (Th) is defined as the number of leaf over the maximum betweenness centrality. The idea behind the Th is that an optimal network should achieve efficient communication while avoiding hub overload. The tree hierarchy has been designed to quantify the balance between both features. Furthermore, we calculated the betweenness centrality (BC) for each region (Freeman, 1977), in order to determine if specific regions differed according to the disability of patients. The BC is defined as the number of shortest paths passing through a given node over the total of the shortest paths of the network (Freeman, 1977). Before moving to the statistical analysis, all the metrics were averaged across epochs in order to obtain one value per subject. An overview of the analysis pipeline is provided in Fig. 1.

Fig. 3

Topological metrics.A. Degree and betweenness centrality: In an undirected, binary network, the degree of a node is defined as the number of links incident upon that node. The higher the degree, the more connected the node in the network. The betweenness centrality (BC) of a node is defined as the number of shortest paths passing through that node over the total number of shortest paths of a network. In the network on the right, it is evident that the red node has a crucial role, although it does not have the highest degree. Such node has the highest BC in the network. B. Leaf fraction and tree hierarchy: a node with degree 1 is defined as a leaf. In the figure, leaf nodes are represented in yellow. On the left, a line-like network, with two leaf nodes. On the right, a star-like network, with a topology that features the maximum possible number of leaf nodes (that is, the total number of nodes minus 1). In line-like networks, it takes on average more steps, as compared to more star-like networks, to go from one node to another. However, the star-like topology relies more heavily on one single node (the central one), that becomes more prone to overload. The tree hierarchy captures the balance between efficient communication and network overload C. Degree correlation: in the upper network, yellow nodes (i.e. lower degree nodes) tend to be linked to red nodes (higher degree nodes). A network where nodes with different degrees tend to form links is defined as disassortative. In the lower network, violet nodes (i.e. lower degree nodes - in this case they are not depicted in yellow as in this example they are not leafs) tend to create links with other violet nodes, while red nodes (with higher degree) tend to link to more red nodes. A network with such topology is defined as assortative. E. Degree divergence: on the x axes there are the degrees of the nodes of a network, on the y axes the frequency with which that value of degree occurred across the network. The broader the range of the degree distribution, the higher the degree divergence.

Statistical analysis

Clinical parameters were compared among the three groups, using ANOVA to test differences in age and educational level and the Chi-square test for gender differences.

For each frequency band, the observations were assumed to be independent. Levene's test was used to test for homogeneity of variances.

The groups were compared for each variable of interest using asymptotic Kruskal-Wallis test in Matlab R2017b (MathWorks®), in order to test the effect of disease stage on brain connectivity in controls, “early stage” and “advanced stage” patients. After the computation of the H-statistic, the p values were corrected using the false discovery rate (FDR) (Benjamini and Hochberg, 1995), so as to account for multiple comparison. The FDR was run across the 90 areas for the BC analysis, and across the four metrics (leaf fraction, degree divergence, degree correlation and tree hierarchy) for the global analysis. For the significant p values (after FDR correction), post-hoc analysis was carried out, using Scheffe correction for multiple comparisons among groups. In all statistical analysis, a level of significance of 0.05 was used.

Furthermore, the disease progression rate (DPR), defined as 48 minus the ALSFRS-R score, divided by the disease duration in months, was used as an index of disease burden. Both the ALSFRS-R and the DPR were related to topological metrics by linear regression. The correlations were corrected across metrics using the FDR.

Results

Table 1 shows the descriptive information of the population. Subject and controls were not different with regard to age (F < 1), education (F < 1) and gender (χ(2) = 1398, p = 0.497). The variables never violated the assumption of equality of variances.

With respect to the tree hierarchy, differences among groups appeared in the delta band [H(2) = 7.66, p = 0.0217, pFDR = 0.0289], in the theta band [H(2) = 10.63, p = 0.0049, pFDR = 0.0065], in the beta band [H(2) = 9.44, p = 0.0089, pFDR = 0.0173] and in the gamma band [H(2) = 10.03, p = 0.0066, pFDR = 0.0266].

In the theta band, post-hoc analysis shows that AS group differed as compared to both ES and HC (p = 0.0082 and p = 0.0476, respectively). The same differences were also present in the beta band ((p = 0.0255 and p = 0.0308, respectively). In the delta band, the controls showed a lower tree hierarchy as compared to AS (p = 0.0262). In the gamma band, a difference has been found between AS and HC as well (p = 0.0112). Results are reported in Fig. 2A.

Fig. 2

Comparison of network parameters by disease stage.A. Scatter plot of the leaf fraction among healthy controls, early stage patients and advanced stage patients in delta, theta, beta and gamma band. B. Scatter plot of the degree divergence among healthy controls, early stage patients and advanced patients in delta, theta, alpha and beta band. C. Scatter plot of the degree correlation among healthy controls, early stage patients and advanced patients in gamma band. D. Scatter plot of the tree hierarchy among healthy controls, early stage patients and advanced patients in delta, theta, beta and gamma band.

The leaf fraction was shown to differ between healthy controls (HC), patients in early stage (ES), and patients in advanced stage (AS) in the delta band [H (2) = 15,88, p < 0.001, pFDR = 0.0014], in the theta band [H (2) = 11,77, p = 0.0028, pFDR = 0.0065], in the beta band [H (2) = 10.02, p = 0.0067, pFDR = 0.0173] and in the gamma band [H (2) = 7.74, p = 0.0208, pFDR = 0.0277]. In the alpha band, the test resulted significant but did not hold to the FDR correction, and it was not tested any further [H(2) = 6.25, p = 0.0438, pFDR = 0.2847].

The post–hoc analysis revealed the same pattern in all four frequency bands, whereby the AS group has higher leaf fraction as compared to ES group and HC in the delta band (p = 0.0235 and p = 0.0005, respectively), in the theta band (p =0.0271 and p=0.0057, respectively), in the beta band (p=0.0149 and p=0.0357, respectively). In the gamma band, only HC and AS were different (p=0.0287). Results are reported in Fig. 2B.

With regard to the degree divergence, differences among groups were evident in the delta band [H (2) = 12,78, p =0.0017, pFDR=0.0034], in the theta band [H (2)=11.35, p=0.0034, pFDR=0.0065], in the alpha band [H(2)=9.72, p=0.0078, pFDR=0.0310], in the beta band [H(2) =8.6856, p=0.0130, pFDR=0.0173] and in the gamma band [H(2)=6.09, p=0.0476, pFDR=0.0476].

The post-hoc analysis revealed that the AS group showed higher degree divergence as compared to both ES and HC groups in the delta band (p = .0367, p = .0026, respectively), theta band (p = .0183, p = .0102, respectively) and beta band (p = .0290 and p = .0490); in the alpha band, AS was only different as compared to HC (p = .0094). In the gamma band, no group was different at the post-hoc analysis, with HC and AS being close to significance (p = .0513). Results are reported in Fig. 2C.

The degree correlation showed differences in the gamma band [H(2) = 8.25, p = 0.0162, pFDR = 0.0277]. Post-hoc analysis revealed that AS differed as compared to HC (p = 0.0199). Results are reported in Fig. 2D.

The BC, computed for the 90 regions of interest, did not show any statistically significant difference in any frequency band after correction for multiple comparisons (data not shown).

The correlation between both the ALSFRS-R and the DPR and the topological metrics did not yield any significant results.

The cumulative distributions of the topological parameters that were different among groups are available in Supplementary material 2.

Discussion

In this paper, we set out to test the hypothesis that the rearrangement of the brain networks in ALS spreads beyond the precentral regions and across multiple frequency bands. Furthermore, we hypothesized that modifications, in terms of increased connectedness of functional brain networks, are related to disease progression.

Applying the PLI and the MST to a well-characterized population of ALS patients compared to HC,

the tree hierarchy did not differ between healthy controls and ES patients. However, both groups differed from AS (see. Fig. 2) in multiple frequency bands. As explained earlier, the tree hierarchy suggest that advanced patients have a suboptimal trade-off between having a well-integrated network while preventing hub overload (see Fig. 3).

Furthermore, the leaf fraction showed a trend, whereby it is the lowest in HC, and then it increases in ES patients and even more so in AS patients. This finding implies a shift toward a more centralized organization of the brain network for patients in the advanced stage as compared to those in the early stage and to controls (Boersma et al., 2013). Similarly, the degree divergence (width of the degree distribution) is lower in controls, and grows in early stages of disease and even higher in patients in more advanced stages. Higher values of degree divergence imply the presence of high-degree nodes (Boersma et al., 2013), which is in line with a more centralized network. Such nodes allow a more rapid synchronization while making the network more vulnerable to targeted attacks (Otte et al., 2015; Tewarie et al., 2015). This happens since, if a high-degree node fails, a high number of nodes of the network become disconnected from each other (Otte et al., 2015).

The provided evidence altogether might be interpreted considering that the less effective trade-off between efficient communication and hub overload (demonstrated by the changes in the tree hierarchy) might be due to a more centralized network that is more prone to node overload (demonstrated by changes in the leaf fraction and degree divergence). Please refer to Fig. 3 for an explanation of the topological metrics.

Furthermore, it has been shown that networks with higher values of degree correlation are more scale–free (Stam and van Straaten, 2012). Recently, it has been proposed that, when a node in the brain network fails, the flow of information that would normally be directed to that node is rerouted to higher nodes. This pattern can be repeated multiple times, propagating across the network and eventually leading to nodes with higher degree to become overloaded (Stam, 2014). Through this process, the network might become more centralized and hyper-connected, maintaining functional output at the cost of becoming more vulnerable to pathological stressors. In accordance with this hypothesis, our results show that the brain functional networks become more star–like (i.e. more connected) and more scale-free as the disease progresses. Furthermore, the analysis of the degree correlation showed that in the gamma band, HC display a less disassortative network, while patients in advanced stage show the highest disassortativity. This implies that when the disease is in an advanced stage, the quantity of connections between nodes with high degree and nodes with low degree is also higher (Newman, 2002). As stated above, this evidence might also be framed within the idea that, as the disease progresses, more and more peripheral nodes need to connect to higher degree nodes, in order to compensate for areas that are no longer able to process information.

Interestingly, the difference in disassortativity was evident in the gamma band. Given that some evidence showed that fast rhythms are spatially restricted (Destexhe et al., 1999; Engel et al., 2013; Steriade et al., 1996), one might speculate that such frequencies would be affected the most by the rewiring process, whereby newer connection across larger (on average) distances would be created in order to connect peripheral areas to more central regions. However, such consideration is only speculative.

Globally, the fact that structurally damaged areas show hyperactivation might be compatible with disease progression relating to more connected networks. One hypothesis might be that hyper connectedness is due to the loss of inhibitory connections (Turner and Kiernan, 2012). We did not find any difference in the centrality of any region (as measured using BC) after correction for multiple comparisons. One fascinating explanation might be that the progression of the pathological process would affect the brain functional networks as a whole more than it does single regions, and hence the global metrics would be more effective in capturing such modifications (Sorrentino et al., 2017). However, given the dependencies among the activity of the brain areas, it could also be possible that the application of the correction for multiple comparisons (FDR) would introduce some false negatives.

Our results were found globally in multiple frequency bands, namely delta, theta, alpha, beta and gamma. This result is in accordance with the previous literature taken altogether. In a MEG study, widespread slow-wave dipole sources could distinguish seven ALS patients from controls. Furthermore, intensified cortical beta desynchronization, followed by delayed rebound, was observed during motor preparation in both contra- and ipsilateral motor cortices of 11 “classic” ALS patients, 9 primary lateral sclerosis patients and 12 asymptomatic carriers of ALS-associated gene mutations compared with healthy control groups (Proudfoot et al., 2017), which is in favour of the hypothesis of cortical hyperexcitability as a key mechanism in the pathogenesis of ALS. Failure of the inhibitory cortical interneuronal function might be underlying such a mechanism (Turner and Kiernan, 2012). Interestingly, Iyer et al. analyzed the EEG of 18 ALS patients and found that connectivity was increased in all frequency bands, and that the clustering coefficient was increased among patients in the alpha and gamma bands (Iyer et al., 2015).

When interpreting our results, one should consider that the choice of using the AAL atlas does affect the results. Classically, the AAL has been considered to have a resolution comparable to that of MEG. Furthermore, it has been reported that the MST itself is affected by the chosen atlas, and applying the MST using the AAL atlas 88% of the individual connections overlapped to the group averaged MST (van Dellen et al., 2018). Possible limitations of our work could be related to patient classification and to the use of a template to reconstruct the sources for subjects lacking MRI scans. With regard to the subsets of patients taken into account, although the classification used may have a low resolution, we still miss better validated biomarkers of disease progression (Turner, 2016). Of note, with regard to the potential usefulness of correlating topological metrics with the currently available clinical scores, we have checked if scores of disability and disease progression, as estimated by, respectively, ALSFRS-R and DPR, were related to MEG topological metrics, revealing no correlation by linear regression analysis after FDR corrections. Regarding to the latter limitation, the quality of the source reconstruction was inferior for subjects' lacking native MRI scans. However, it was recently shown that for the metrics we chose, no statistically significant difference was evident in the same subject reconstructed with either native or template MRI (Douw et al., 2017).

Conclusions

In conclusion, our study is the first magnetoencephalography study on a large cohort of ALS patients. We estimated the connectivity by measuring synchronization between brain areas (instead of coactivation) and used a graph theoretical approach to study disease staging in ALS. The value of our study is to have shown that the stage of the disease in ALS is related to a widespread topological reorganization of the brain toward a more integrated and more vulnerable network.