Disrupted core-periphery structure of multimodal brain networks in Alzheimer’s disease

In Alzheimer’s disease (AD), the progressive atrophy leads to aberrant network reconfigurations both at structural and functional levels. In such network reorganization, the core and peripheral nodes appear to be crucial for the prediction of clinical outcome because of their ability to influence large-scale functional integration. However, the role of the different types of brain connectivity in such prediction still remains unclear. Using a multiplex network approach we integrated information from DWI, fMRI, and MEG brain connectivity to extract an enriched description of the core-periphery structure in a group of AD patients and age-matched controls. Globally, the regional coreness—that is, the probability of a region to be in the multiplex core—significantly decreased in AD patients as result of a random disconnection process initiated by the neurodegeneration. Locally, the most impacted areas were in the core of the network—including temporal, parietal, and occipital areas—while we reported compensatory increments for the peripheral regions in the sensorimotor system. Furthermore, these network changes significantly predicted the cognitive and memory impairment of patients. Taken together these results indicate that a more accurate description of neurodegenerative diseases can be obtained from the multimodal integration of neuroimaging-derived network data.


INTRODUCTION
The brain is a complex network where differently specialized areas are anatomically and func-Complex network: A collection of elements that interact dyadically in a nontrivial way. tionally connected. Because of such interconnected structure, focal damages can affect the rest of the network through the interruption of communication pathways. Indeed, many neurological disorders affecting language, motor, and sensory abilities are often due to a disconnection syndrome caused by the anatomical connectivity breakdown between the relevant brain areas (Geschwind, 1965;Schmahmann & Pandya, 2008). In the case of neurodegenerative diseases, the disconnection hypothesis is supported by a progressive death of neurons and synapses that induce gross atrophy. Empirical evidence has shown that Alzheimer's disease (AD) patients with severe motor and cognitive impairments exhibited anatomical disconnections among regions between cerebral hemispheres that resemble those observed in split-brain subjects (Delbeuck, Collette, & Van der Linden, 2007;Lakmache, Lassonde, Gauthier, Frigon, & Lepore, 1998). In Parkinson's disease (PD) intrahemispheric dissociations between subcortical and cortical structures have been linked to disturbances in cognition, perception, emotion, and sleep (Cronin-Golomb, 2010). In addition, functional connectivity alterations within and between hemispheres have been reported in both AD (Adler, Brassen, & Jajcevic, 2003;Babiloni et al., 2009;Blinowska et al., 2016;Sankari, 2010) and PD (Dubbelink et al., 2013;Luo et al., 2015), suggesting their potential role in early diagnosis.
Altogether, these findings suggest that neurodegenerative diseases should be considered as a network problem. Recent approaches based on network theory have greatly advanced our understanding of the connection mechanisms characterizing brain diseases (Stam, 2014). Among others, decreased efficiency, modularity, and hub centrality have been largely reported in neurodegeneration and associated with the stage of disease. Increasing evidence suggests that the core-periphery structure of the human connectome-supporting global integration Connectome: A network description of a neural system's wiring; an exhaustive list of the physical connections (e.g., synapses, projections, fiber tracts) that link all neural elements (e.g., neurons, neuronal populations, macroscopic brain regions). of information among distant areas-is highly affected by the AD process and that resulting changes might be effective predictors of cognitive declines. On one hand, brain areas forming the core of the network-that is, central and mutually connected nodes-have been reported to be preferentially attacked by AD (Yan et al., 2018). On the other hand, brain regions forming the periphery of the network-that is, nodes that are only weakly connected to the other units in the network-appear to be crucial for the degeneration (Daianu et al., 2015). While these results refer to structural brain connectivity, the relative contribution of functional brain connectivity into the network core-periphery changes remains poorly understood.
Based on the aforementioned empirical and theoretical grounds, we hypothesize that neurodegeneration would affect the core-periphery structure of the brain network at both anatomical and functional levels. More specifically, we expected that the extraction of the core-periphery organization by integrating information from multimodal brain networks would give more accurate predictors of AD and cognitive impairment. Finally, based on the evidence that hubs are the most attacked nodes in many neurological diseases and psychiatric disorders (van den Heuvel & Sporns, 2013), we hypothesize that the core brain regions would be mostly impacted by the AD atrophy process.
To test these predictions, we considered multiple brain networks derived from DWI, fMRI, and MEG data recorded in a group of AD patients and age-matched healthy controls (HC; Figure 1. Multiplex brain network construction. Different neuroimaging data are collected and preprocessed separately. We used the Desikan cortical atlas parcellation (Desikan et al., 2006) to infer connectivity networks from DWI, fMRI, and MEG source-reconstructed data. The color of the line indicates the software that has been used in each step of the pipeline. We spatially aligned all the estimated brain networks to construct the multiplex brain network. Figure 1). Cognitive impairments in AD patients were described using multidomain behavioral measurements. We extracted the multimodal core-periphery structure of the brain networks through a multiplex network approach, where all the available information is kept at differ-Multiplex network: A multilayer network where interlayer connections are allowed only between homologous nodes. ent connectivity layers. Multiplex network theory has been recently introduced to specifically model and analyze complex systems whose units can be linked through different types of connectivity (Battiston, Nicosia, & Latora, 2014;Boccaletti et al., 2014;De Domenico et al., 2013;Kivelä et al., 2014). Main applications in neuroscience have focused on the characterization of higher order network motifs (Battiston, Nicosia, Chavez, & Latora, 2017;Crofts, Forrester, & O'Dea, 2016;Pedersen, Zalesky, Omidvarnia, & Jackson, 2018) and node centrality (Bentley et al., 2016;De Domenico, Sasai, & Arenas, 2016;Guillon et al., 2017) in both healthy and diseased conditions (De Domenico et al., 2016;Guillon et al., 2017;Yu et al., 2017). Here, we evaluated how AD impacted the multiplex core-periphery organization (Battiston, Guillon, Chavez, Latora, & De Vico Fallani, 2018), and we tested the correlation of the regional coreness with the cognitive and memory impairment of patients. See the Material and Methods section for more details on the experimental design and methods of analysis.

Multimodal Core of Brain Networks
We integrated multimodal information by constructing nine-layer multiplex brain networks containing DWI, fMRI, and MEG connectivity between 68 cortical regions of interest (ROIs; Material and Methods). To estimate the likelihood of each ROI i to be in the multiplex core we computed its coreness C i by counting how many times it was in the multiplex core across different density thresholds (Battiston et al., 2018). At each threshold, the multiplex core-periphery structure was obtained by linearly combining the node strength of all the layers through a vector parameter c (Material and Methods).
Because we do not know a priori the best combination, we derived the optimal c * by using a data-driven approach that efficiently explores the parameter space to maximize the difference between AD and HC regional coreness. Specifically, we used the particles swarm optimization algorithm (PSO) to maximize the Fisher's criterion F(c) (Material and Methods). Results show that the optimal c * components are found to be highly heterogenous and that the DWI layer, as well as MEG-alpha1 and fMRI layers, are the main contributors to separate the AD and HC group (Table 1, Figure 2A).
In the HC group, the multiplex core tended to include large portions of temporal, superior parietal, and occipital cortices, and to a minor extent central and superior frontal regions ( Figure 2B). On average AD patients exhibited a loss of coreness with respect to HC particularly Figure 2. Regional coreness of the multiplex brain networks. Panel A), shows the results of the particle swarm optimization (PSO) used to find the best layer coefficients vector c that maximizes the Fisher score F(c) between AD and HC subjects. In the upper plot, each dot represents the position of a particle at a given iteration in the original 9-dimension coreness contribution coefficient vector space. The color of the dots code for the corresponding Fisher score. Results were projected over the three main network layers for the sake of illustration The other nonshown components were rapidly zeroing-out during the 81 iterations needed to converge to the optimum as shown in the bottom plot. Panel B) shows the corresponding average coreness for the healthy control (HC) population and for the Alzheimer's disease (AD) group. The position of the nodes identifies the barycenter of each ROI in the cortical surface, here represented in transparent gray; the color of each node codes for the average corenessC i . Figure 3. Differences in regional coreness between AD and HC subjects. Panel A) shows the between-group difference of corenessC AD,i − C HC,i as a function of the healthy population's corenessC HC,i ; the slope of the regression line in gray measures the coreness disruption index κ = −0.20. The color of the circles code for the difference between average coreness in the AD and HC groups; stars point out the ROIs for which we reported a significant difference (p < 0.025, Supplementary Table S2; Guillon et al., 2019). Panel B) illustrates the values of the between-group coreness difference over the Desikan cortical atlas. Color code is the same as in panel A.
in the temporal, superior parietal, and occipital cortices. These regions were already known to form the core of multiplex brain networks derived from DTI and fMRI data (Battiston et al., 2018).

Reorganization of Core-Periphery Structure in AD
To quantify the observed network changes, we defined the coreness disruption index κ as the slope of the line obtained by regressing the difference between the average coreness (at each ROI and across subjects) of the two groups with the average coreness of the healthy one (Termenon, Achard, Jaillard, & Delon-Martin, 2016) (Material and Methods). We found a significant negative κ value, indicating that AD preferentially attack ROIs with a high coreness (κ = −0.20, p = 2.45e−10). This result was also consistent at the individual level when we extracted the coreness disruption index in each patient (Supplementary Table S1; Guillon et al., 2019). In particular, by statistically comparing the average coreness of the two groups, we reported a significant decrease of coreness in core regions, such as temporal, parietal, and occipital cortices as well as a significant increase of coreness in the right paracentral area that are instead more peripheral (p < 0.025, Figure 3A, B, Supplementary Table S2; Guillon et al., 2019).
Based on the hypothesis that AD is a disconnection syndrome (Delbeuck et al., 2007;Geschwind, 1965) leading to disorganized network configurations (Sanz-Arigita et al., 2010), we next generated a series of synthetic multiplex networks starting from the ones observed in the HC group and attacking an increasing amount of links in each layer. Specifically, we simulated the dysconnection process by decreasing the weights of the links that were preferentially connected to the core nodes (Materials and Methods). Results show that the coreness disruption index decreased with the number of links that were randomly attacked when the reduction was sufficiently strong, that is, 75%. In this situation, we could generate the same κ values observed in the multiplex brain networks of the AD group by attacking between 25% and 30% of the links ( Figure 4A). Notably, this result could not be obtained when we attacked the links preferentially connected to peripheral nodes ( Figure 4B). Altogether, these findings indicate that AD is associated with a pervasive random reconfiguration of brain connectivity that primarily affects the nodes of the multiplex core.

Coreness Disruption Predicts Cognitive and Memory Deficits
We finally conducted a correlation analysis to better understand how the observed multiplex brain network changes were associated with the behavioral performance of AD patients. Results show that both cognitive and memory deficits could be predicted by the individual loss of regional coreness. At the global scale, the coreness disruption index significantly correlated with the Mini-Mental State Examination (MMSE) (R = 0.46, p = 0.028) as well as with the immediate (R = 0.47, p = 0.024) and free recall (R = 0.59, p = 0.005) scores. The higher the κ values, the better was the performance of the patients ( Figure 5A, Supplementary material; Guillon et al., 2019). At the local scale, temporal, parietal, and cortices were highly positively correlated with the behavior of patients. Notably, these ROIs overlapped with those exhibiting significant decreases of regional coreness with respect to healthy controls ( Figure 3B). We found similar positive correlations for bilateral middle frontal ROIs (R = 0.36, p = 0.092 for left, R = 0.35, p = 0.100 for right), while areas in the motor system appeared not to be involved except for the paracentral lobule that tended to negatively correlate with the MMSE (R = −0.55, p = 0.007) and immediate recall scores (R = −0.36, p = 0.089; Figure 5B).

Multiplex Brain Networks
The increasing availability of multimodal neuroimaging data holds a great potential to enrich our knowledge about fundemental neural mechanisms and to improve the precision of predictive biomarkers of brain diseases (Calhoun & Sui, 2016). However, how to integrate information from different neuroimaging modalities is still an open issue. Existing approaches have mainly focused on merging information at the level of the native data structure (e.g., signal or images; Biessmann, Plis, Meinecke, Eichele, & Muller, 2011;Uludag & Roebroeck, 2014). Only recently, investigators have started to propose fusion algorithms in an effort to infer brain connectivity (Ng, Varoquaux, Poline, & Thirion, 2012) or to detect mental states (Lei et al., 2011). Here, we adopted a complementary solution-based on the nascent field of multilayer network theory-that preserves the original nature of the different connectivity types. Sim-Multilayer network: A collection of elements that interact dyadically according to different types of connectivity. Each connectivity type is represented in a different layer, so that both intra-and interlayer connections are allowed. ilar approaches have been already used in the case of temporal (Bassett et al., 2011;Pedersen et al., 2018), multifrequency (Brookes et al., 2016;De Domenico et al., 2016;Guillon et al., 2017;Tewarie et al., 2016;Yu et al., 2017), and DTI-fMRI brain networks (Battiston et al., 2017) in human but also in nonhuman species (Bentley et al., 2016;Crofts et al., 2016). This study considers for the first time brain networks obtained from three different neuroimaging modalities-DWI, fMRI, and MEG-to construct multiplex brain networks consisting of nine connectivity layers and to derive an augmented description of their core-periphery structure in healthy and Alzheimer's diseased subjects. A crucial step in the characterization of multiplex networks is how to weight the contribution from different layers (Boccaletti et al., 2014;Kivelä et al., 2014), which typically contain connectivity measured in different units (e.g., number of fiber tracks and amount of signal correlation). While this is in general an arbitrary choice, here we established an objective way to associate a weight to each layer by maximizing the difference of regional coreness-that is, the likelihood of each region to be in the core-between the groups. Results showed that all three modalities are necessary to the group separation. In particular, for MEG only alpha1 was determinant while the other frequency layers had very low, or null, weights. This is in line with current evidence showing that the alpha1 frequency band contains the most discriminant power and connectivity changes in AD (Babiloni et al., 2004;Blinowska et al., 2016). Notably, DWI had a very high contribution coefficient as compared with the other layers. Core-periphery structure of diffusion-based networks is known to be very robust (Hagmann et al., 2008;van den Heuvel & Sporns, 2011) with respect to functional layers, and this might possibly depend on the heterogeneity of the weighted node degree distribution. More in general, further research is needed to derive criteria that objectively balance the layer contribution in the extraction of multiplex network properties.

Network Reorganization in Alzheimer's Disease
AD is associated with network changes affecting the structure and function of the brain at multiple spatial and temporal scales (Stam, 2014). It has been hypothesized that these network reconfigurations could result from dysconnection patterns initiated by the gross atrophy of the brain. While several studies have found significant changes in terms of network efficiency, modularity, and node centrality, the direction of these alterations-in terms of increments or decrements with respect to healthy controls-is often unclear and modality dependent (Tijms et al., 2013). Here, we focused on the core-periphery structure of the human brain, which has been shown to have a significant impact on cognition ensuring global integration across remote cortical areas (van den Heuvel & Sporns, 2011). Structural connectome studies have reported that AD patients, from the preclinical to dementia stages, have significant hub-concentrated lesion distributions (Brier et al., 2014;Buckner et al., 2009;Crossley et al., 2014;Dai et al., 2014;Shu, Wang, Bi, Zhao, & Han, 2018). However, recent evidence is suggesting that network disruption is prevalent in the peripheral network components in both AD (Daianu et al., 2015) and mild cognitive impairement (MCI) patients (Zhao et al., 2017). These inconsistent findings suggest that the network disruption mechanisms remain unclear. By integrating information from structural and functional brain networks, we aimed to overcome this controversy and provide a more comprehensive insight. Our multiplex network approach shows that core regions were globally affected in AD patients as compared with HC subjects and that this result could be modeled by a global random rewiring process. Specifically, we reported significant decrements of coreness in temporal and parietal cortices, which are heavily affected by atrophy processes and beta-amyloid deposition (Buckner, Andrews-Hanna, & Schacter, 2008). However, this change was paralleled by a significant increase of coreness in the paracentral lobules, which originally belonged to the multiplex periphery. Because regions of the sensorimotor system-such as paracentral lobule-are not directly affected by the atrophy process (Agosta et al., 2010), we speculate that possible compensatory mechanisms could have therefore taken place. In line with this hypothesis, recent findings suggest that more efficient motor commands in mild cognitive impaired patients could trigger the later functional decline (Kubicki, Fautrelle, Bourrelier, Rouaud, & Mourey, 2016). Longitudinal studies involving healthy subjects converting into AD will be fundamental to confirm or reject this prediction (Dubois et al., 2016).

Connectivity-Based Biomarkers of Clinical Behavior
Brain wiring organization is critically associated with human cognition and behavior as well as with several neurological and psychiatric disorders (Stam, 2014). Network indices describing core-periphery and rich-club organization in structural brain networks have been shown to predict cognitive and motor deficits in multiple sclerosis (Stellmann et al., 2017), and Huntington disease (Harrington et al., 2015), as well as communication impairment in schizophrenia (van den Heuvel et al., 2013). More pertinent to this work, rich-club biomarkers extracted from DTI networks have been shown to correlate with cognitive and memory deficits in Alzheimer's disease (Daianu et al., 2015;Stam et al., 2009;Tijms et al., 2013).
Here, we showed that the coreness disruption index-quantifying the global tendency to weaken core-periphery structure in multimodal brain networks-determined the cognitive and memory performance of our AD patients. Patients with a stronger core-periphery organization had higher MMSE and Free and Cued Selective Reminding Test (FCSRT) scores. At the local scale, temporal, parietal, as well as frontal areas tended to positively correlate with patients' behavior. These association regions have been shown to be implicated in the prediction of AD cognitive performance (Khachiyants & Kim, 2012) and more in general in memory and language (Gordon, 1995;Pochon et al., 2002;Squire, 1992). We also found negative correlations with the paracentral lobule (especially right), a region that is typically involved in motor-related tasks but not in integrative functions.
From a network perspective, the coreness of regions that tended to be in the multiplex core-such as temporal, parietal, and occipital cortices-were positively correlated with patients' performance, while among the peripheral areas the paracentral lobule was negatively correlated with the behavior. This means that in the presence of more severe cognitive and memory deficits, the relative decrease of connectivity in core regions tended to be replaced by periphery components of the brain system. This result would confirm the existence of an adjusting mechanism, where the sensorimotor system might be involved in the compensation of connectivity loss in systems that are directly impacted by amyloid-beta plaques and tau neurofibrillary tangles accumulation (Iaccarino et al., 2018).

Methodological Considerations
The basic algorithm behind the detection of the core-periphery structure in multiplex networks is purely deterministic (Ma & Mondragón, 2015). This means that in principle we could not evaluate the statistical relevance of the identified structure. To overcome this limitation, we adopted a procedure that consisted in extracting the core-periphery structure from a series of multiplex networks obtained by filtering the actual brain multiplex network with increasing density thresholds (Battiston et al., 2018). This way we could derive a probabilistic measure of coreness by counting how many times each ROI was assigned to the core across all the possible thresholds. For the sake of simplicity we filtered each brain multiplex by retaining the strongest links so that the average node degree of each layer ranged from k = 1 to k = N − 1 (Material and Methods).
After filtering we did not binarize the surviving links so that we applied the core-periphery algorithm to sparse weighted multiplex networks. This approach allows us to exploit all the available information in the multiplex brain networks. At the same time, we remark that additional care is needed, as it introduces issues related to the different nature and distributions of the link weights (Buldú & Papo, 2018). Here, we mitigated this problem by using in the core-periphery algorithm the vector c of parameters that can weight the contribution of each layer (Material and Methods). Alternative solutions have been recently proposed taking into account the normalization of the weight across the layers by means of singular value decomposition (Mandke et al., 2018). Finally, because of the difficulty to measure connectivity between different neuroimaging modalities we did not consider interlayer links in our multiplex representation, as instead explored in more theoretical studies (Buldú & Porter, 2017;De Domenico et al., 2016). Despite the absence of interlinks, our multiplex network approach is able to extract higher order topological properties that cannot be obtained by other single-layer approaches. Notably, given the high level of overlap and correlations between the layers, our multiplex networks cannot be reduced to simple colored-edge graphs, where different colors Graph: A mathematical description of a network; elements and their interactions are represented as nodes/vertices and connections/edges, respectively. are associated with different types of connectivity (Boccaletti et al., 2014).
We used an optimization algorithm-namely the particle swarm optimization-to find the best combination of c components that maximized the difference of the coreness between AD and HC subjects (Material and Methods). This method presents two limitations that are important to mention here. First, the time complexity increases exponentially with the number of layers M in order to find a stable solution. We verified that for M > 10 the research complexity becomes rapidly intractable because of the large space of parameter combination to explore. Second, the cost function optimized by the algorithm and used to evaluate how segregated the two groups are (i.e., two sets of coreness vectors) should be carefully chosen as its accuracy is highly impacted by the size of the feature space (here N, the number of ROIs) and the size of the samples (here the size of the cohorts). More advanced techniques taking into account the possible nonlinear and/or non-Euclidean nature of the feature space should be considered for very large networks (e.g., support vector machines, Riemannian geometry).

Conclusion
Consistent with our hypothesis, we have shown that the AD atrophy process generates multimodal connectivity changes that can be quantified by a multilayer network approach. Specifically we have identified that both core and-to a minor extent-peripheral cortical areas are affected in AD, and that the direction of the effect was opposite. Decrease of coreness in temporal, parietal, and occipital areas-forming the rich core of the human brain-is paralleled by a possible compensatory increment in cortical regions that are in the sensorimotor system and that are more peripheral. These cortical network signatures varied over individuals and were significant predictors of cognitive and memory deficits. Furthermore, we reported a general framework for the statistical comparison of core-periphery organization in arbitrary multiplex networks. Taken together, our results offer new insights into the crucial role of core-periphery organization in neurodegenerative diseases.

Cohort Inclusion
The study involved 23 Alzheimer's diseased (AD) patients (13 women) and 26 healthy agematched control (HC) subjects (19 women). All participants underwent the Mini-Mental State Examination (MMSE) for global cognition and the Free and Cued Selective Reminding Test (FCSRT) for verbal episodic memory. Inclusion criteria for all participants were (a) age between 50 and 90; (b) absence of general evolutive pathology; (c) no previous history of psychiatric diseases; (d) no contraindication to MRI examination; and (e) French as a mother tongue. Specific criteria for AD patients were (a) clinical diagnosis of Alzheimer's disease; and (b) Mini-Mental State Examination (MMSE) score greater or equal to 18. All subjects gave written informed consent for participation in the study, which was approved by the local ethics committee of the Pitie-Salpetriere Hospital. All experiments were performed in accordance with relevant guidelines and regulation.

Data Acquisition and Preprocessing
Magnetic resonance imaging (MRI) acquisitions were obtained using a 3T system (Siemens Trio, 32-channel system, with a 12-channel head coil). The MRI examination included (a) 3D T1-weighted volumetric magnetization-prepared rapid gradient echo (MPRAGE) sequence with the following parameters: thickness = 1 mm isotropic, repetition time (TR) = 2,300 ms, echo time (TE) = 4.18 ms, inversion time (TI) = 900 ms, acquisition matrix = 256 256; (b) echo planar imaging (EPI) sequence with the following parameters: one image with no diffusion sensitization (b0 image) and 50 diffusion-weighted images (DWI) at b = 1,500 s/mm 2 , thickness = 2 mm isotropic, TR = 13,000 ms, TE = 92 ms, flip angle = 90, acquisition matrix = 128 116; (c) functional MRI (fMRI) resting-state sequence sensitive to blood oxygenation leveldependent (BOLD) contrast with the following parameters: 200 images, thickness = 3 mm isotropic, TR = 2,400 ms, TE = 30 ms, flip angle = 90, acquisition matrix = 64 64. All MR images were processed using the Clinica software (http://www.clinica.run). We first used the t1-freesurfer-cross-sectional pipeline to process T1-weighted images. This pipeline is a wrapper of different tools of the FreeSurfer software (http://surfer.nmr.mgh.harvard.edu/; Fischl, 2012). It includes segmentation of subcortical structures, extraction of cortical surfaces, cortical thickness estimation, spatial normalization onto the FreeSurfer surface template (FsAverage), and parcellation of cortical regions. Functional MRI images preprocessing has been conducted using the fmri-preprocessing pipeline. Slice timing correction, head motion correction, and unwarping have been applied using SPM12 tools (www.fil.ion.ucl.ac. uk/spm). Separately, the brain mask has been extracted from the T1 image of each subject using FreeSurfer. The resulting fMRI images have then been registered to the brain-masked T1 image of each subject using SPM's registration tool. Finally, diffusion-weighted images have been processed using the dwi-preprocessing pipeline of Clinica. For each subject, all raw DWI volumes were rigidly registered (6 df ) to the reference b0 image (DWI volume with no diffusion sensitization) to correct for head motion. The diffusion-weighting directions were appropriately updated (Leemans & Jones, 2009). An affine registration (12 df ) was then performed between each DWI volume and the reference b0 to correct for eddy current distortions. These registrations were done using the FSL flirt tool (www.fmrib.ox.ac.uk/fsl). To correct for EPI-induced susceptibility artifacts, the field map image was used as proposed by Jezzard and Balaban (1995) with the FSL prelude/fugue tools. Finally, the DWI volumes were corrected for nonuniform intensity using the ANTs N4 bias correction algorithm (Tustison et al., 2010). A single multiplicative bias field from the reference b0 image was estimated, as suggested in Jeurissen, Tournier, Dhollander, Connelly, and Sijbers (2014).
The magnetoencephalography (MEG) experimental protocol consisted in a resting state with eyes closed (EC). Subjects were seated comfortably in a dimly lit electromagnetically and acoustically shielded room and were asked to relax. MEG signals were collected using a whole-head MEG system with 102 magnetometers and 204 planar gradiometers (Elekta Neuromag TRIUX MEG system) at a sampling rate of 1,000 Hz and online low-pass filtered at 330 Hz. The ground electrode was located on the right shoulder blade. An electrocardiogram (EKG, Ag/AgCl electrodes) was placed on the left abdomen for artifacts correction and a vertical electrooculogram (EOG) was simultaneously recorded. Four small coils were attached to the participant in order to monitor head position and to provide coregistration with the anatomical MRI. The physical landmarks (the nasion, the left and right preauricular points) were digitized using a Polhemus Fastrak digitizer (Polhemus, Colchester, VT). We extracted three consecutive clean epochs of approximately 2 min each.
Signal space separation was performed using MaxFilter to remove external noise. We used in-house software to remove cardiac and ocular blink artifacts from MEG signals by means of principal component analysis. We visually inspected the preprocessed MEG signals in order to remove epochs that still presented spurious contamination. At the end of the process, we obtained a coherent dataset consisting of three clean preprocessed epochs per subject. We reconstructed the MEG activity on the cortical surface by using a source imaging technique (Baillet et al., 2001;He, 1999): (a) We used the previously segmented T1-weighted images of each single subject (Fischl et al., 2002(Fischl et al., , 2004 to import cortical surfaces in the Brainstorm software (Tadel et al., 2011) where they were modeled with approximately 20,000 equivalent current dipoles (i.e., the vertices of the cortical meshes). (b) We applied the wMNE (weighted minimum norm estimate) algorithm with overlapping spheres (Lin et al., 2006) to solve the linear inverse problem. Both magnetometer and gradiometer, whose position has been registered on the T1 image using the digitized head points, were used to localize the activity over the cortical surface.

Construction of Brain Networks
We built, for each modality, one or multiple brain connectivity networks whose nodes are regions of interest (ROIs) defined by the Desikan cortical atlas parcellation (Desikan et al., 2006; N = 68 regions); and links are weighted by a given connectivity measure estimated between each pair of nodes resulting in 68 × 68 fully symmetric adjacency matrices.
In the case of MEG, we used the spectral coherence as a connectivity estimator with the following parameters: window length = 2 s, window type = sliding Hanning, overlap = 25%, number of FFT points (NFFT) = 2,000 for a frequency resolution of 0.5 Hz between 2 Hz and 45 Hz included.
For fMRI data, we focused our analysis on the scale 2 wavelet correlation matrices that represented-with a TR = 2,400ms-the functional connectivity in the frequency interval 0.05-0.10Hz (Achard, 2006;Bassett & Bullmore, 2009;Biswal, Yetkin, Haughton, & Hyde, 1995;Cordes et al., 2001;De Vico Fallani, Richiardi, Chavez, & Achard, 2014). This choice was mainly motivated by the fact that the interpretation of different frequencies in fMRI is not clearly defined (Chen & Glover, 2015), whereas in E/MEG specific mental states can be more directly associated with distinct bands.
For DWI data, we used the Clinica software to estimate the fiber orientation distributions (FODs) using constrained spherical deconvolution (CSD) algorithm from MRtrix3 dwi2fod tool and tractography based on iFOD2 algorithm from MRtrix3 tckgen tool. The connectome is finally estimated by counting the number of tracts connecting each pair of nodes according to the given parcellation file using MRtrix3 tck2connectome tool.

Network Methods and Models
We constructed multiplex brain networks in each subject by spatially aligning DWI, fMRI, and MEG source reconstructed connectivity networks. This led to the following multiplex network with M = 9 layers: are respectively the largest and smallest entries of W [m] . This way, all the links' weights ranged between 0 and 1 and became comparable quantities across layers.
To extract the coreness of the nodes from the resulting multiplex networks, we followed the procedure described by Battiston et al. (2018). First, we filtered each layer by preserving the strongest weights for a broad range of increasing thresholds. Specifically, we considered density-based thresholds so that each layer had the same average node degree from k = 1 to k = N − 1. Then, for each threshold we computed the core-periphery of the filtered multiplex network by evaluating (a) the multiplex richness μ i of node i, defined as follows: with s [m] i the strength of the node in the m-th layer, and c [m] the components of the vector c that modulate the contribution of each modality-specific layer. And (b), similarly to the original paper, we decomposed the richness function into two components based on the links of node i that are going towards nodes with lower richness and those towards nodes with higher richness s [m] = s [m]− + s [m]+ . Thus, the multiplex richness of a node towards richer nodes is defined as follows: We finally counted the number of times that each node was in the core across all the explored thresholds, and we normalized by the maximum theoretical value. As a result, we obtained the coreness C i that can be written as follows: where δ [k] i = 1, if node i is in the core for the average node degree k. 0, otherwise.
To simulate the disconnection process, we generated random multiplex networks by decreasing the intensity of the links in each layer starting from the actual multiplex brain network of the HC group. First, we fixed the percentage of links to be attacked and the percentage of weight reduction. Then, we randomly attacked the links by selecting those with higher probability to be connected to core or periphery nodes. For each HC individual, we generated n rand new randomized multiplexes according to the following pseudocode: Step 1. Initialization (a) Fix the number of links to randomly attack (L) (b) Fix the percentage of weight reduction (R) (c) Give a probability p(i) to each node proportional to its coreness Step 2. Repeat until termination criteria are met For k = 1, ..., L do (a) Pick a random node i with a probability p(i) (b) Pick a random node j with a probability p(j) (c) Decrease the weight of the link w ij by R (d) Exclude the link i − j in the next iterations (e) k ← k + 1 Step 3. Normalize and output connectivity matrix We chose the minimum number of randomizations necessary to obtain a variance approximately equal to the one observed in the HC and AD groups. This number was n rand = 3 and gave in total N RA = n rand × N HC = 78 samples.

Particles Swarm Optimization and Statistical Analysis
We used the PSO algorithm (Kennedy & Eberhart, 1995) under the MATLAB(R) software with the default parameters. The Fisher's criterion F(c) was defined as follows: withĪ Pop (c), the average local (i.e., node level) index, here the coreness C, over a population Pop, which in our case belongs to {AD, HC}, and, s 2 with s a subject belonging to the population Pop.
Since, in our case, F(c) = F(ac), ∀a ∈ R + , and in order to save one dimension in the searching space, we expressed the coefficient c as a point on the positive section of the unitary hypersphere of dimension M = 9 such that sin φ 1 ... sin φ 8 sin φ 1 ... sin φ 7 cos φ 8 sin φ 1 ... sin φ 6 cos φ 7 sin φ 1 ... sin φ 5 cos φ 6 sin φ 1 ... sin φ 4 cos φ 5 sin φ 1 ... sin φ 3 cos φ 4 sin φ 1 sin φ 2 cos φ 3 sin φ 1 cos φ 2 To consider the non-Gaussian nature of the data we considered nonparametric statistics when assessing differences between populations and prediction of behavioral scores. To these ends, we used respectively permutation t tests and Spearman correlation coefficients. The statistical threholds were set to α = 0.05, and we applied a rough false discovery rate (FDR) correction to account for the N = 68 post hoc tests at the level of brain regions (α FDR = 0.025).