## PET-MR data

67 healthy volunteers (mean age: 51.1 ± 16.4 years, range 20–82 years, almost uniformly distributed) were recruited prospectively between December 2015 and February 2017. The main exclusion criteria consisted of major internal pathology, diabetes mellitus, cancer, absence of a first-degree relative with dementia, history of important neurological and/or psychiatric disorders or substance abuse or pre-study use of centrally acting medication. All subjects underwent a complete neurological examination, performed by a board-certified physician, had a mini-mental state examination (MMSE) score ≥ 28 and their index on the Becks Depression Inventory (BDI) was ≤ 9. The study was approved by the Ethics committee of the University Hospital Leuven (study number S58571 – Belgian Registration Number B322201526273) and was conducted in full accordance with the latest version of the Declaration of Helsinki. All participants provided written informed consent before inclusion in the study.

18F-FDG was administered by intravenous injection of 151.9 ± 9.8 MBq. All subjects underwent simultaneous FDG PET and MR scanning on a hybrid 3T Signa PET-MR scanner (GE healthcare, Chicago, IL, USA). List mode images were acquired upon tracer injection in the scanner for 60 minutes, from which static (40–60 min pi) data were derived. The first 15 minutes of the simultaneous scan, no MR sequences were applied in order not to invoke primary auditory cortex activation. MR image acquisition was performed using an 8-channel phased-array coil. In addition to a whole brain volumetric T1-weighted image (3D BRAVO, TR/TE = 8.5/3.2 ms, 1x1x1 mm voxel size) and a fluid-attenuated inversion recovery (FLAIR) image (3D CUBE, TR/TE = 8500/130 ms, 1x1x1.4 mm voxel size) were collected. Resting-state data in eyes-open condition were acquired with TR/TE = 1.7s/2ms, flip angle = 90, voxel size = 3.6 x 3.6 x 4 mm. DTI and reverse phased DTI were acquired using a b-value of 1500 s/mm², applied along 48 uniformly distributed directions (TR/TE = 12000/85 ms, 2.5 mm isotropic voxel size mm). Default vendor-based MRAC (MR-based attenuation correction, v1.0) corrected PET images were reconstructed using ordered subset expectation maximization (OSEM) with 6 iterations and 28 subsets, and post-smoothed with a 3 mm isotropic Gaussian filter.

Patients showing too much movement were excluded from the start, since this can have a large impact on rs-fMRI data and subsequent analysis, with exclusion criterium: mean framewise displacement higher than 0.3 mm.

To investigate the effect of ageing, we selected a “young” (\(n=26\)) and “elderly” (\(n=28\)) population by setting an age threshold at age ≤ 45 and age ≥ 55 respectively and identified the core nodes for both groups.

## Brain connectivity

Both structural, functional and metabolic networks consisted of \(100\) cortical nodes, defined by regions of interest (ROIs) obtained with a Schaefer parcellation scheme [16]. As such, each network was represented by a \(100\times 100\) adjacency matrix, describing the connectivity or edge weights between each pair of nodes. SC and FC connectomes for a specific population were calculated as the average weighted network across the population. For each modality, the different type of information and different way of network construction resulted in an unequal level of network sparsity of the corresponding weighted adjacency matrix. Therefore, all networks were binarized to ensure that the connectivity density, i.e. the total number of edges, of each network was equal. However, since the choice of binarizing threshold is rather arbitrary, results were averaged across a full range of connectivity densities (from 10% up to 50%,stepsize: 1%).

**Structural connectivity**

Diffusion images were processed using MRtrix3 [17] and the FMRIB Software Library [18]. Preprocessing of the diffusion MRI data included denoising, Gibbs ringing removal, correction for EPI susceptibility, eddy-current-induced distortions, gradient-nonlinearities and subject motion. From the corrected diffusion data, response functions for single-fiber white matter (WM), grey matter (GM) and cerebrospinal fluid (CSF) were estimated using an unsupervised method [19]. Single-shell 3-tissue constrained spherical deconvolution (SS3T-CSD) was performed to obtain the fiber orientation distributions (FODs) for WM, GM and CSF [20] using MRtrix3Tissue (https://3Tissue.github.io), a fork of MRtrix3 [17]. Consequently, the FODs were used to conduct probabilistic tractography using the Fiber Assignment by Continuous Tracking (FACT) algorithm, which tracks the trajectory of white matter tracts by propagating streamlines along the primary direction of water diffusion at each voxel [21]. Anatomically Constrained Tractography was performed alongside FACT using the tissue-segmented T1-weighted image in order to ensure that the generated streamlines were biologically accurate [22]. Whole-brain tractograms were re-weighted using Spherically Informed Filtering of Tractograms 2 (SIFT2) [23], which adjusted the streamline weights in order to represent the underlying diffusion signal more accurately. For each subject, the tractogram and Schaefer parcellation were combined to produce a subject-specific structural connectome. The corresponding edge weights were defined by the number of streamlines between two nodes, normalized by the volumes of both regions represented by the two nodes such that each value of the SC matrix reflected the density of the white matter streamlines between the corresponding two nodes.

**Functional connectivity**

After pre-processing the BOLD rs-fMRI time-series for each individual, edge weights for the functional networks were calculated by the Pearson correlation coefficient between the average time-series in each ROI. Preprocessing of the BOLD time-series was performed using fmriprep version 1.5.10 [24]. Each T1-weighted (T1w) volume was corrected for intensity non-uniformity using N4BiasFieldCorrection [217] and skull-stripped using Advanced Normalization Tools (ANTs) [25]. Spatial normalization to the ICBM 152 Nonlinear Asymmetrical template (version 2009c) [26] was performed through nonlinear registration with the ANTs registration tool using brain-extracted versions of both T1w volume and template. Brain-tissue segmentation of CSF, WM and GM was performed on the brain-extracted T1w volume using FSL [27]. Functional data were slice-time corrected using AFNI (3dTshift) [28], and realigned to a mean reference image using FSL (mcflirt) [29]. Fieldmap-less distortion correction was performed by co-registering the functional image to the intensity-inverted T1w image [30], constrained with an EPI distortion atlas [31] and implemented with ANTs (antsRegistration). This was followed by co-registration to the corresponding T1w volume using boundary-based image registration [32] with nine degrees of freedom. Framewise displacement was calculated for each functional run [33]. The non-aggressive variant of ICA-based Automatic Removal Of Motion Artifacts (AROMA) was used to generate and remove noise components from the fmriprep-processed output [34]. Subsequently, mean WM, mean CSF and global mean signals were calculated and regressed out in a single step using least squared regression. Finally, the data were filtered with a high-frequency bandpass filter of [0.01Hz, 0.1Hz] to exclude confounding high-frequency content.

**Metabolic connectivity**

For each subject, the 18F-FDG uptake was normalized by proportional scaling, i.e. dividing the uptake by the total uptake in the grey matter, after smoothing the data with a Gaussian kernel with a Full Width Half Maximum (FWHM) of 8mm. Although different strategies exist to define metabolic connectivity networks, all these approaches rely on the assessment of regional co-variation in 18F-FDG uptake across subjects, which is different from the subject-specific approaches to define structural and functional connectivity. We applied Sparse Inverse Covariance Estimation (SICE) [35], also known as Gaussian graphical modelling. Basically, SICE imposes a sparsity constraint on the maximum likelihood solution of the inverse covariance (IC) matrix under the assumption of a Gaussian model, which means that the sample size can be less than the number of brain regions modeled. Since the brain network organization is assumed to be sparse [36], SICE is considered as a valid approach to model metabolic brain connectivity. Although SICE has proven to be an effective tool for identifying the structure of an IC matrix, its use is not recommended to estimate the magnitude of the non-zero entries in case of weighted adjacency matrices [35]. However, since we work with binary instead of weighted matrices, SICE is an appropriate choice for our approach to estimate the zero and non-zero entries of the IC matrix.

## Multilayer network

A multilayer network consists of several classical networks, each layer encoding a specific type of information. Since multilayer networks can no longer be represented by classical adjacency matrices, we used a tensor formalism to describe these networks. More specifically, we defined a multilayer adjacency tensor of \(N\) nodes and \(L\) layers with components given by \({M}_{i\alpha }^{j\beta }\), encoding the connectivity between node \(i\) in layer \(\alpha\) and node \(j\) in layer \(\beta\) (\(i,j=\text{1,2},\dots ,N\); \(\alpha ,\beta =\text{1,2},\dots ,L\)). To easily represent the tensor, we used the standard approach of flattening this rank-4 tensor into a rank-2 tensor, also known as the supra-adjacency matrix, with the diagonal blocks encoding the inter-layer connectivity. For our three-modal PET-MR data, the multilayer network consisted of three layers, i.e. a structural, functional and metabolic connectivity layer, resulting in a SC – FC – MC multilayer network. In order to naturally extend classical network metrics, the different layers should be interconnected since otherwise the analysis of the multilayer adjacency tensor comes down to analyzing each layer separately. To combine individual layers, links are added between corresponding nodes across layers, either ordinally by linking corresponding nodes between adjacent layers only, or categorically, by linking a node to the corresponding nodes across all layers. We implemented the latter option since ordinally linking assumes that neighboring networks are prioritized.

Since core nodes tend to score high on nearly all centrality measures, we focused on a combination of degree centrality (DC) and eigenvector centrality (EC) for defining the core nodes in the multilayer networks. In a single layer network with adjacency matrix \(A\), the DC of node \(i\) is given by its number of connections, where the EC of node\(i\) represented by \({x}_{i}\) is defined as:

\({x}_{i}=\frac{1}{\lambda }{\sum }_{k}^{}{A}_{ki}{x}_{k}\) (6.1)

As a result, EC measures a node’s importance while taking into account the importance of its neighbors. This can be written in matrix form \(Ax=\lambda x\), which is the eigenvector equation of \(A\). Since the adjacency matrix \(A\) is positive definite, the Perron-Frobenius theorem states that there exist a unique and positive eigenvector corresponding to the leading eigenvalue \({\lambda }_{1}\), yielding the EC values for all nodes. EC values were calculated by the power iteration method, which makes sense intuitively. All nodes start with equal EC, but as the computation progresses, nodes with more edges start gaining importance, which mainly propagates to the nodes to which they are connected. After several iterations, EC values stabilize, resulting in the final EC estimates.

Both DC and EC measures for nodes in monoplex network were generalized towards multilayer networks. The (scalar) degree of a node in a multilayer network is either given by the degree of the aggregated topological network, i.e. the union of the monoplexes, or by the overlapping (summed) degree \({o}_{i}={\sum }_{\alpha }^{}{k}_{i}^{\alpha }\), where \({k}_{i}^{\alpha }\) is the degree of node \(i\)in layer \(\alpha\). Since both degree measures tend to be highly correlated, the overlapping degree centrality (ODC) was used. On the other hand, EC generalization is less trivial and there are several ways to do so [37]. The most elegant way is to rewrite the eigenvalue equation using an equivalent tensor formulation (single layer):

\({A}_{i}^{j}{x}_{j}={\lambda }_{1}{x}_{i}\) (6.2)

using the Einstein summation convention \({A}_{i}^{j}{x}_{j}\equiv {\sum }_{j}^{}{A}_{i}^{j}{x}_{j}\). As a result, the adjacency matrix is now given by the rank-2 tensor \({A}_{i}^{j}\) containing one contravariant and one covariant component. The generalization of the above equation to the multilayer case with adjacency tensor \({M}_{i\alpha }^{j\beta }\) is then given by:

\({M}_{i\alpha }^{j\beta }{\varTheta }_{j\beta }={\lambda }_{1}{\varTheta }_{i\alpha }\) (6.3)

with \({\lambda }_{1}\) the leading eigenvalue and \({\varTheta }_{i\alpha }\) the corresponding eigentensor. This equation can easily be solved using the supra-adjacency matrix formulation:

\(\left(\begin{array}{ccc}{M}^{1}& w{I}_{N}& w{I}_{N}\\ w{I}_{N}& {M}^{2}& w{I}_{N}\\ {wI}_{N}& w{I}_{N}& {M}^{3}\end{array}\right)\left(\begin{array}{c}{{\Theta }}_{}^{1 }\\ {{\Theta }}_{}^{2 }\\ {{\Theta }}_{}^{3 }\end{array}\right)={\lambda }_{1}\left(\begin{array}{c}{{\Theta }}_{}^{1 }\\ {{\Theta }}_{}^{2 }\\ {{\Theta }}_{}^{3 }\end{array}\right)\) (6.4)

with \({I}_{N}\) the unity matrix of size \(N\times N\),\({M}^{\alpha }\) corresponds to the single-layer adjacency matrix of layer \(\alpha\), \({\varTheta }_{i}^{\alpha }\equiv {\varTheta }_{i\alpha }\)encodes the eigentensor centrality (ETC) of each node \(i\) in layer \(\alpha\) while accounting for the whole interconnected structure, and \(w\) is an intra-layer weight factor. The ETC value \({\theta }_{i}\) of each node is found by contracting \({\varTheta }_{i}\)with the rank-1 tensors \({u}^{\alpha }\) with all components equal to 1, i.e. \({\theta }_{i}={\varTheta }_{i\alpha }{u}^{\alpha }\). The choice of this aggregation corresponds to a maximum entropy principle, which is a valid choice when all layers are considered equally important [38]. The weight factor \(w\) is chosen such that the total number of inter-layer connections is equal to the intra-layer connections, i.e.

\(w=\frac{3{\gamma C}_{N}^{2}}{6N}=\frac{\gamma (N-1)}{4}\) (6.5)

where \(N\) is the number of nodes, \({C}_{N}^{2}\) is the number of 2-combinations out of N, i.e. \({C}_{N}^{2}=N(N-1)/2\), and \(\gamma\) is the connectivity density of the monoplex networks.

## Core selection

At each binarizing threshold between 0.1 an 0.5 (step 0.01), we identified a set of core nodes by selecting these nodes scoring high, i.e. a fraction \(\delta\) of the standard deviation \(\sigma\)above the mean value \(\mu\) (\(\mu +\delta *\sigma\)), on both DC/ODC and EC/ETC (depending on the monoplex/multiplex nature of the network). Since the definition of “scoring high” depends on the parameter \(\delta\), and therefore is rather arbitrary, we calculated a set of core nodes for different \(\delta\) values (from 0.4 to 1.6, stepsize 0.2). In this way, the coreness coefficient \({C}_{i}\) of each node \(i\) is calculated as the normalized frequency (across binarizing thresholds and \(\delta\) values) of being part of the core. To make the interpretation more convenient, we colored the top 15% of the most central nodes in a different color.

To quantify the similarity between the cores of layer \(\alpha\) and layer \(\beta\), the core similarity coefficient is defined as:

\({S}_{c}=\frac{{I}_{c}^{\alpha \beta }}{{N}_{c}^{\alpha }} (0< {S}_{c}<1)\) (6.6)

where \({I}_{c}^{\alpha \beta }\) is the number of nodes which are part of the core of both layer \(\alpha\) and \(\beta\), and \({N}_{c}^{\alpha }\) is the total number of nodes in layer \(\alpha\) [39].

First, the structural (SC), functional (FC) and metabolic (MC) connectomes are calculated for the average population. Consequently, the corresponding monoplex cores are identified based on the combined (DC, EC) measure, and the multiplex core is obtained by selecting nodes scoring high on both ODC and ETC. SC and FC population connectomes were obtained by calculating the average weighted SC and FC networks respectively.