Networks of myelin covariance

Abstract Networks of anatomical covariance have been widely used to study connectivity patterns in both normal and pathological brains based on the concurrent changes of morphometric measures (i.e., cortical thickness) between brain structures across subjects (Evans, 2013). However, the existence of networks of microstructural changes within brain tissue has been largely unexplored so far. In this article, we studied in vivo the concurrent myelination processes among brain anatomical structures that gathered together emerge to form nonrandom networks. We name these “networks of myelin covariance” (Myelin‐Nets). The Myelin‐Nets were built from quantitative Magnetization Transfer data—an in‐vivo magnetic resonance imaging (MRI) marker of myelin content. The synchronicity of the variations in myelin content between anatomical regions was measured by computing the Pearson's correlation coefficient. We were especially interested in elucidating the effect of age on the topological organization of the Myelin‐Nets. We therefore selected two age groups: Young‐Age (20–31 years old) and Old‐Age (60–71 years old) and a pool of participants from 48 to 87 years old for a Myelin‐Nets aging trajectory study. We found that the topological organization of the Myelin‐Nets is strongly shaped by aging processes. The global myelin correlation strength, between homologous regions and locally in different brain lobes, showed a significant dependence on age. Interestingly, we also showed that the aging process modulates the resilience of the Myelin‐Nets to damage of principal network structures. In summary, this work sheds light on the organizational principles driving myelination and myelin degeneration in brain gray matter and how such patterns are modulated by aging.

subject of profound changes across the lifespan with gray matter myelination following inverted-U shaped trajectories with age (Grydeland, Walhovd, Tamnes, Westlye, & Fjell, 2013) whose characteristics exhibit specific anatomical distribution patterns (Dean et al., 2016). So far, however, the important question of the interactions between local cortical myelination changes remains largely unexplored. Do myelin changes occur independently across brain regions? Is myelination in one region related to or modulated by variations in other regions? Are synchronized myelination changes across the cortical mantle topologically organized? Could common underlying factors (e.g., genetics) be driving these concurrent myelin fluctuations? Our paper is precisely motivated by such questions.
The analysis of the covariance between regional estimates of morphometric features derived from magnetic resonance images (MRI) is well established (Andrews, Halpern, & Purves, 1997;Mechelli, Friston, Frackowiak, & Price, 2005). Recent studies have investigated in detail the topological organization of the brain (He, Chen, & Evans, 2007;Sanabria-Diaz et al., 2010) to provide evidence for a network of anatomical covariance in the healthy and diseased brain (Alexander-Bloch, Giedd, & Bullmore, 2013;Evans, 2013). A number of studies have shown strong age-related effects on the number and strength of the interhemispheric correlations in cortical thickness and three key network topological properties: characteristic path length, clustering index, and local efficiency (Chen, He, Rosa-Neto, Gong, & Evans, 2011;Wu et al., 2012;Yang, Tsai, Liu, Huang, & Lin, 2016;Zhu et al., 2012).
Morphological covariance studies are based on macroscopic features of the brain and do not provide an insight into the concurrent changes taking place at the microscopic scale within brain tissue. The latter changes are of primary interest in neuroscience but are largely intractable from the most widely used types of anatomical MRI scans (e.g., T1-weighted, T2-weighted). A recent study has shown that microscopic changes within brain tissue could lead to the spurious detection of apparent morphological change (Lorio et al., 2016). New methods have emerged based on the ratio of T1-weighted and T2-weighted images (Glasser and Van Essen, 2011) that provide an insight into the biological processes underlying brain tissue changes. However multiple microscopic properties of the brain may impact ratio estimates, complicating the interpretability of observed findings. Quantitative MRI (qMRI) addresses the limitations of standard MRI anatomical data by providing quantitative estimates of the MRI parameters that drive signal intensities in an MRI image. The correlation between these estimates and microscopic features of brain tissue has motivated their use as in vivo markers of microstructure (Lutti, Dick, Sereno, & Weiskopf, 2014). qMRI data is corrected for the sources of artifact that affect standard anatomical MRI images, leading to enhanced reproducibility and sensitivity to physiological brain changes (Weiskopf et al., 2013).
Because myelin and iron concentrations are the main contributors to qMRI estimates (Fukunaga et al., 2010;Schmierer et al., 2007;St€ uber et al., 2014), the qMRI changes reported in neuroscience applications have mainly been attributed to iron deposition and demyelination processes, in line with histological findings (Callaghan et al., 2014;Draganski et al., 2011). qMRI also allows the in vivo delineation of the heavily myelinated boundaries of visual (Sereno, Lutti, Weiskopf, & Dick, 2013) and primary auditory (Dick et al., 2012) areas.
This study is based on the MRI parameter, Magnetization Transfer (MT) (Helms, Dathe, Kallenberg, & Dechent, 2008b), which exhibits a high level of specificity toward tissue myelination (Callaghan et al., 2014;Helms, Draganski, Frackowiak, Ashburner, & Weiskopf, 2009;Lorio et al., 2014;Lorio et al., 2016). The MT mechanism is based on the exchange of magnetization between free water and protons bound to macromolecules. MT values are mainly driven by the local macromolecule density and the amount of water in close proximity with these macromolecules. Amongst the macromolecules involved, myelin cholesterol has been suggested as a major contributor to magnetization transfer (Koenig, Brown, Spiller, & Lundbom, 1990;Koenig, 1991).
Recently, using MT maps, Hunt et al. (2016) evidenced the presence of myelin concurrent changes among brain anatomical structures but limited to studying its predictive value for electrophysiological functional connectivity.
In this study, we use in vivo measures of MT to investigate patterns of correlations in myelination change between gray matter regions. In particular, we aim to study how myelination, a key microstructural feature of brain tissue, changes synchronously among spatially distant regions to form characteristic networks of myelin covariance (Myelin-Nets). Based on the graph theory framework we explore, for the first time, the topological organization of the Myelin-Nets and its modulation by age. Two age groups (Young-Age, 20-31 years old; Old-Age, 60-71 years old) were compared to describe agerelated Myelin-Nets topological changes in an elapsed time of 40 years.
The continuous aging trajectory of the topology of the Myelin-Nets was also investigated after the fifth decade of life. This approach could be considered as another step for revealing basic principles of gray matter organization, and how these are modified by aging processes.

| Participants
The dataset included 562 participants of the CoLaus/PsyCoLaus cohort (Firmann et al., 2008;Preisig et al., 2009) and other research studies carried out at the LREN laboratory. The age range was 18-87 years (277 Females). In our first study, a pool of 151 participants were selected to define two age groups: Young-Age adults: 73 (35 Females), comprising all participants of the dataset from 20 to 31 years old (mean age 5 24.56 years) and Old-Age adults: 78 (40 Females), from 60 to 71 years old (mean age 5 65.04 years) (see Figure 3, Results).
The age groups were selected with a gap of 40 years to capture different lifespan stages of the myelination processes.
Owing to the reduced number of participants in our cohort between 30 and 47 years old, our second study, aiming to explore the continuous aging trajectories of the Myelin-Nets, was limited to the 48-87 years old age range, which involved a subset of 437 out of the 562 subjects of our original data.

| Ethics statement
The study participants gave written informed consent at the time of their enrollment and completed questionnaires approved by the local MELIE-GARCIA ET AL.

| 1533
Ethics Committee. The authors state that they have obtained approval from CHUV and Colaus Data Sharing and Publications Committee for use of the data and confirm that the data was analyzed anonymously.

| MRI acquisition and preprocessing
Participants were examined on a 3 T whole-body MRI system (Magnetom Prisma, Siemens Medical Systems, Germany), using a 64-channel RF receive head coil and body coil for transmission. On visual inspection study participants showed neither macroscopic brain abnormalities, that is, major atrophy, nor signs of overt vascular pathology (i.e., microbleeds and white matter lesions. Participants with extended atrophy or with white matter hyperintensities (WMH) of grade 2 or more by the Scheltens rating scale (Scheltens et al., 1993) were not included.
The whole-brain quantitative protocol comprised three multi-echo 3D fast low angle shot (FLASH) acquisitions with predominantly Magnetization Transfer-weighted (MTw: TR/a 5 24.5 ms/68), proton density-weighted (PDw: TR/a 5 24.5 ms/68) and T1-weighted (24.5 ms/218) contrast (Helms et al., , 2008a(Helms et al., , 2008b. The MTw contrast was achieved by use of a Gaussian-shaped RF pulse prior to the excitation (4 ms duration, 2208 nominal flip angle, 2 kHz frequency offset from water resonance). Multiple gradient echoes were acquired with alternating readout polarity with a minimal echo time TE of 2.34 ms and a time span of 2.34 ms between consecutive echoes; 6/8/8 echoes were acquired for the MTw/PDw/T1w acquisitions to keep the TR value identical for all acquisitions. The image resolution was 1 mm 3 isotropic, the field of view was 256 3 240 3 176 mm and the matrix size was 256 3 240 3 176. Parallel imaging was used along the phase-encoding direction (acceleration factor 2 GRAPPA reconstruction (Griswold et al., 2002)), 6/8 partial Fourier was used in the partition direction. The acquisition time was 7 min per contrast. Data were acquired to calculate maps of the radio frequency (RF) transmit field B11 using a 3-D echoplanar imaging (EPI) spin-echo (SE) and stimulated echo (STE) method (Lutti, Hutton, Finsterbusch, Helms, & Weiskopf, 2010;Lutti et al., 2012) and to correct for effects of RF transmit inhomogeneities on the quantitative maps Helms and Dechent, 2009;Weiskopf et al., 2013). The image resolution of the B1-mapping data was 4 mm 3 .
The echo time was 39.06 ms and TR was set to 500 ms. The spin-echo flip angle was decreased from 2308 to 1308 in steps of 108 . The acquisition time of the B1 mapping data was 3 min. B0-field mapping data was acquired using a 2-D double-echo FLASH sequence to correct for geometric distortions in the 3-D EPI data as described in Lutti et al. (2010Lutti et al. ( , 2012 ucl.ac.uk/spm) and Matlab 7.11 (Mathworks, Sherborn, MA, USA). The MT maps were computed as described in Helms, Dathe, and Dechent (2008a) and Helms and Dechent (2009) using the MTw, PDw, and T1w images with minimal echo time (TE 5 2.34 ms) to minimize R2* bias on the MT estimates (Lorio et al., 2016). The amplitude of the MT effect is governed by the duration, power and off-resonance frequency of the MT saturation pulse. Therefore, the MT effect is most often characterized using semi-quantitative measures, that is, that depend on the set of acquisition parameters kept constant for all study participants. Note that these MT measures differ from the common Magnetization Transfer Ratio (MTR), by accounting for local T1 relaxation and flip angle inhomogenity effects, resulting in enhanced robustness and sensitivity to myelin concentration (Helms et al., 2008b). In short, the Neuromorphometrics atlasing methodology (Figure 1, step 2) consists of two main steps. First, each individual MT image is segmented into three different brain tissue classes (cerebral spinal fluid, gray matter and white matter) using the "Segment" (unified segmentation) tool in SPM12, which includes a registration to MNI (Montreal Neurological Institute) space. Second, the probabilistic atlas of each of the 114 anatomical structures of the Neuromorphometrics atlas (see further) are spatially registered with the extracted gray and white matter tissue maps using the "Shoot" tool in SPM12, based on a nonlinear advanced registration algorithm (Ashburner and Friston, 2011). Rules of probability are used to properly combine the previous images to ultimately obtain a probabilistic label map for each brain structure.

| Construction of the myelin data matrix
At every gray matter voxel (in subject space), the probability of belonging to a specific anatomical structure is provided. From above, a maximum probability label maps are calculated at all gray matter voxels (in subject space) which are labeled according to the structure of maxi-

| Construction of the myelin correlation matrix
The process of constructing the myelin correlation matrix is summarized in Figure 2. For each age group, a linear regression was performed on the regional MT estimates to remove the effects of age, age 2 , Therefore, the global MT was not regressed out as confounding variable. Global MT effect subtraction in principle may highlight interregional differences but at the same time will affect genuine covariance patterns (and topological network attributes) without the possibility of estimating the introduced bias (Borchardt et al., 2016;Fox, Zhang, Snyder, & Raichle, 2009;Murphy, Birn, Handwerker, Jones, & Bandettini, 2009;Murphy and Fox, 2017;Schwarz and McGonigle, 2011).
We defined a connection as the statistical association in MT values (a surrogate measure of myelin) between each pair of brain regions for a parcellation scheme of 114 anatomical structures. The statistical similarity or synchronized co-variations in myelination between two regions was measured by computing the Pearson's correlation coefficient across subjects (Figure 2, steps c and d) (He et al., 2007). Hence, the interregional Pearson's correlation matrix C (N 3 N, N (5114) being the number of gray matter regions)-or Myelin correlation matrix-contains elements C ij that are the value of the Pearson's correlation between regions i and j. Self-connections were excluded, implying zeros in the diagonal of this symmetric matrix.
It is important to point out that a partial correlation analysis could not be used in our case because the sample size was not large enough for a robust estimation of this measure (i.e., the number of structures in the Neuromorphometrics parcellation is higher than the number of subjects for each group).
In the next step, we obtained for each age group Nboot 5 2000 bootstrap samples of the connectivity matrix by selecting a random subset of subjects with replacement using the classical bootstrapping procedure described in (Efron and Tibshirani, 1994). This algorithm guarantees independent samples. Bootstrapping allows the computation of empirical distributions of the connectivity matrices and of the topological network properties (i.e., mean, standard deviation, etc.) that, with a sufficient number of data points, reflect the true underlying distributions. This could not be obtained with a permutation resampling approach. In addition, and very, important the bootstrapped connectivity matrices allowed us to estimate the significance of changes in network properties between age groups (Young-Age and Old-Age) taking into account the variability of these properties within each group. The connectivity matrices obtained from the 2000 bootstraps were thresholded to create sparse binary graphs (Figure 2, step e). Rather than restricting our analysis to a binarized graph obtained by applying a single threshold value, we explored the Network Properties of the graphs over a range of sparsity degrees varied from 0.5 to 0.9 in steps of 0.02 (a sparsity degree of 0.9 means that 90% of the connectivity matrix is discarded, keeping only the highest 10% of the connectivity values) (Figure 2 Simmons, 2011;Sanabria-Diaz et al., 2010). This range of sparsity degree was chosen to allow for all network properties to be accurately estimated and the number of spurious edges in each network minimized as indicated in previous studies (Achard and Bullmore, 2007;He et al., 2007).
Corresponding threshold values R k were calculated for each of the 2000 correlation matrices so that their elements C ij were set to 1 when |C ij | > R k and 0 otherwise. This procedure normalizes the networks to have the same number of nodes and edges, enabling the examination of the relative network properties obtained for each group.
2.6 | Graph analysis to characterize the "networks of myelin covariance" (Myelin-Nets) A great number of natural systems can be represented by complex networks. Graph Theory is usually considered an attractive model for the mathematical treatment of such networks, including those representing brain connectivity (Sporns, 2011). In general, a complex network can be represented as a graph G 5 [N,K], where the nodes N are the components of this system and the edges (K), the relations or connections between them (Boccaletti, Latora, Moreno, Chavez, & Hwang, 2006). In our specific case, the nodes are the anatomical regions obtained from the brain parcellation and the edges are the co-variations (correlation) in myelination across subjects between pairs of these brain regions.
It is important to note here that this is a mathematically derived network, whose connections do not explicitly reflect anatomical or physiological mechanisms in the brain. However, because these "Networks of myelin covariance" (Myelin-Nets) are based on anatomical data, they can be taken as possible biomarkers of underlying biophysical mechanisms. These Myelin-Nets are unweighted because all edges are assumed to indicate relations of equivalent strength between nodes, and undirected, simply summarizing symmetric relations (such as correlations) between nodes.
We used graph theory to study the myelin co-variation networks in the Young-Age and Old-Age groups. This mathematical treatment allows us to characterize the age-related changes of global and local phenomena observed when myelination in any structure fluctuates concurrently with myelination in its neighborhoods and other distant brain regions of the network. In other words, graph theory gives us a framework to explore the Myelin-Nets architecture and how efficiently the information of myelin fluctuations is "exchanged" over the network (in terms of the graph theory). Importantly, these networks cannot be interpreted in terms of temporal causality. First, because myelination co-variations are assessed across subjects and second, because we are using Pearson correlations, measures that do not provide directional/ causal information on the interactions between pairs of nodes.
In particular, we analyzed the following global network attributes: cluster index, local and global efficiency, and characteristic path length.
To describe the nodal properties of the network we computed the betweenness centrality attribute that allowed us to identify the network hubs. Additionally we carried out a "Targeted Attack" study to evaluate the resilience of the Myelin-Nets when the most important regions (hubs) are virtually attacked. In the following, these measures will be defined with the traditional interpretation of general networks.
However, their usefulness as relevant descriptors of anatomical brain states will become apparent in the next sections.

| Clustering index (C)
The clustering index C i of a node "i" is defined as the number of existing connections between the node's neighbors divided by all their possible connections. It is a measure of the inherent tendency of nodes to cluster into strictly connected neighborhoods (Watts and Strogatz, 1998). Nodes are considered neighbors when a connection between them exists, which is not reduced to a physical neighborhood concept.
The clustering index for the whole graph G is defined as the average clustering around each node: represents the number of nodes. Clearly, 0 < C < 1; and C 5 1 if and only if the network is fully connected; that is, each node is connected to all other nodes.

| Characteristic path length (L)
The characteristic path length L of the graph G is the smallest number of connections required to connect one node to another, averaged over all pairs of nodes. It is a measure of the typical separation between two nodes (structures) i and j 8i; j 2 N ð Þ , and it is defined as the mean of geodesic lengths d ij over all pairs of nodes.
In the unweighted network context, the geodesic length d ij is defined as the number of edges along the shortest path connecting nodes i and j (Boccaletti et al., 2006;Watts, 1999;Watts and Strogatz, 1998).

| Network efficiency
The concept of efficiency has also been expressed in terms of information flow (Latora and Marchiori, 2001). That is, small world networks are very efficient in terms of global and local communication and they are defined to have high global E glob and local E loc efficiency. The global E glob of a graph G is expressed as: This measure reflects how efficiently the information can be exchanged over the network, considering a parallel system in which each node sends information concurrently along the network. On the other hand, the E glob of G is defined as the average efficiency of the local subgraphs: where G i is the subgraph of the neighbors of "i." This measure reveals how much the system is fault tolerant, showing how efficient the communication is among the first neighbors of i when it is removed (Latora and Marchiori, 2001). As above, nodes are considered neighbors when a connection between them exists, which is not reduced to a physical neighborhood concept.
2.10 | Nodal centrality: Normalized betweenness centrality (NBC) The "betweenness centrality" B i of a node i is defined as the number of shortest paths between any two nodes that run through node i (Freeman, 1977). We measured the normalized betweenness centrality as

| Methodology for studying differences in the myelination correlation across brain lobes
To investigate differences in myelination correlation (interconnectivity) between age groups across brain lobes we used the anatomical subdivision of the brain in lobes proposed by Tzourio-Mazoyer et al. (2002). We assessed the intra lobe myelination connectivity as the mean of the abso-  (He et al., 2009;Wu et al., 2011b). It is worth noting that the monotonic changes of the NP curves with the sparsity degree make the area under the curve a suitable descriptor of the global performance of the networks.
We followed three main steps to examine differences in network properties between groups: (a) construction of the empirical boot- percentile bootstrap CI) (Efron, 1982)

| Age modulation of the Myelin-Nets global network properties: influence of the gray matter parcellation
There is empirical evidence that topological properties of brain networks depend on the gray matter parcellation used (Zalesky et al., 2010). In order to study the effect of grey matter parcellation on the

| Methodology for studying the aging trajectory of the Myelin-Net's global network attributes
To uncover dynamical properties of the Myelin-Nets topological organization with age, we proposed a continuous aging trajectory analysis based on a sliding window approach. This methodology has been used in two forms to study networks of anatomical covariance: (a) using a weighted contribution of the data points in the windows to the correlation coefficient at the age window centers (Zalesky et al., 2015); (b) using a fixed number of data points (subjects) per window, the overlap between contiguous windows being selected heuristically (Vasa et al., 2017). Here, we kept the number of subjects per window constant, equal to the number of subjects of the Young-Age group (73), to keep the correlation matrices and their topological features across windows unbiased by the number of data points. For each slid, the youngest subject of the current window was replaced with the nearest older participant (i.e., step size 5 1: two contiguous windows only differed by one subject). This process was repeated iteratively across the age range 5 48-75 years old (Section 2.1). For each sliding-window, the Myelin-Nets and their topological attributes were calculated using the methodology described in Construction of the Myelin correlation matrix section without bootstrapping. The sliding window "age" was defined as the median age of the participants in each window. The minimum "age" difference between contiguous windows was 0.01 (4 days) and the maximum was 0.6 (half a year). The maximum age difference between the youngest and oldest subjects within a window was 9.16 years (around the "window age" 5 74 years old), the minimum 3.29 years and the mean 5.3 years.
We focused our attention on the main global network attributes of the Myelin-Nets: clustering index, characteristic path length, local and global efficiency, global connectivity, and the connectivity strength between homologous regions. Additionally a study of the aging trajectory of the normalized betweenness centrality (NBC), as nodal network property, was performed.
The aging trajectory of the Myelin-Nets attributes were fitted as a function of the "window age" using a polynomial model. The order of the polynomial fitting was determined by the Akaike's Information Criterion (AIC) (Bozdogan, 1987). The statistical significance of the polynomial coefficients was assessed through Student's t test in the linear regression model.

| RE S U L TS
3.1 | Global gray matter myelination changes with age Before detrending for age, visualization of the full 562 dataset showed that the global myelination in gray matter followed an inverted-U shape trajectory with age ( Figure 3).
This motivated the use of a second-order polynomial model to fit the age dependence of the MT data for the detrending (using "polyfit. m" subroutine implemented in MATLAB 2015a), with significant coefficients (p < .05). Aging was found to have a heterogeneous effect on myelination across gray matter regions (Supporting Information, Figures   S1-S4). The peak in myelination was found to take place between 40 and 50 years of age for most of brain anatomical regions (Supporting Information, Figure S2). The modeled age-dependence of the regional MT values was detrended from the original values before calculation of the correlation matrices.
3.2 | Aging modulates the correlation strength between homologous regions but not the global myelin correlation strength

| Aging modulates the correlations in myelination in brain lobes
As shown in Figure 5, the correlation in myelination within brain lobes (pulled from both hemispheres) is differently modulated by age. The Old-Age group showed higher intra-lobe correlation strength than the Young-Age group except for the temporal, parietal, and occipital lobes, where no significant differences were found (the detailed statistical results can be found in Supporting Information, Table SII).

| Most connected structures in Myelin-Nets
In a more detailed study of the correlations in myelination, we found a set of structures where the variations in myelination were particularly strongly correlated with myelination changes in other gray matter regions. Figure 6 shows the 15 structures with the strongest mean inter-regional covariations in each age group (the cortical surface plots were created using the BrainNet Viewer package (http://www.nitrc. org/projects/bnv) (Xia, Wang, & He, 2013)). Table 1

| Effects of age on the Myelin-Nets hubs
We also studied the effects of age on the hubs of the Myelin-Nets, defined as having a Normalized Betweenness Centrality (NBC, a quantitative measure of the importance of structures in the Myelin-Nets) over 1.5 (Figure 7, panel A).
The full list of hub structures can be found in Table 2 In contrast, the Young-Age group depicted higher global efficiency ( Figure 9, panel a). Details of the results of the statistical tests can be found in Supporting Information, Table SIV. 3.7 | Effects of age on the resilience of the myelin networks to targeted attack In the "Targeted Attack" study, we found that the resilience to virtual damage of the principal structures in the Myelin-Nets increases with age: the Old-Age group showed higher resilience after "simulated attacks" to the Myelin-Net hubs ( Figure 10).
In contrast, the Young-Age group was more vulnerable than the    Table 3 shows the summary of the network property (AUC) statistics and targeted attack study, comparing both age groups, for the the parcellations. Interestingly, the AAL atlas was also reported to yield the lowest homogeneity in resting state activity in Gordon's study (Gordon et al., 2016). In contrast, the Neuromorphometrics atlas provided the highest sensitivity to differences between age groups in all  NPs except for the global connectivity attribute. It is important to point out that clustering index, characteristic path length, and local and global efficiency were the most reproducible measures across brain parcellations. Another relevant and consistent result is that global connectivity showed no significant differences between age groups independently of the brain parcellation.

| D I SCUSSION
In this study, we provide empirical evidence that myelination processes in gray matter are spatially correlated, irrespective of the distance between brain regions. When pooling together all possible covariations between pairs of structures, a nonrandom network of concurrent changes emerged: the network of myelination covariance (Myelin-Net).
The Myelin-Net showed well-structured topological patterns strongly  The structures common to both groups are highlighted in bold font.

| Correlated myelination processes in the brain: Possible mechanisms
In this article, we found that myelination fluctuates concurrently across the brain gray matter with well-defined patterns. Several mechanisms have been proposed to explain the concurrent fluctuations between the spatially distinct regions of the networks of anatomical covariance: (a) common normal/pathological vulnerabilities; (b) neurodevelopment; (c) genetic factors (Alexander-Bloch et al., 2013); and (d) the presence of axonal connections between brain structures (Gong, He, Chen, & Evans, 2012). We hypothesize that these putative mechanisms may explain the synchronicity of myelination processes between gray matter regions.
One could expect that if two brain structures are anatomically connected by nervous fibers, and if the density of myelin along the fiber trajectories is relatively constant and distributed across axonal branches inside each region, these distant structures might show correlated myelin densities. If a proportion of nervous fibers were to degenerate, myelin density in the gray matter structures connected by these fibers might concurrently decrease. Alternatively, intracortical myelin changes may act to optimize synchronicity and timing between structures and their reciprocal connections (Haroutunian et al., 2014). These homeostatic mechanisms may be responsible for correlated myelination processes, acting to preserve an optimal conduction delay (timing) and thus synchrony of neural information processing between connected structures (Kimura and Itami, 2009;Salami, Itami, Tsumoto, & Kimura, 2003). Hunt et al. (2016) and Huntenburg et al. (2017)   observation and attempts to explain how the brain reacts to different challenges. This theory represents the brain as a dynamically adaptive structure that changes in both positive and negative ways with age.
During neurodevelopment and normal aging, myelination processes are expected to correlate among groups of structures as a consequence of the establishment and consolidation of neural circuits in the former case, and to compensate aging-related neural degradations (i.e., cortical thinning, myelin deterioration, etc.) in the latter. It is known that the degradation of myelin (presence of balloons, etc.) across the lifespan triggers continuous remyelination as repair mechanism (Peters, 2002).
This phenomenon can be thought of as one that would increase correlations in myelination between structures if they are involved in the same pathological events. Also the need to recruit additional neural resource to sustain cognitive performance in middle-age and old-age brains increases the chance of new myelination correlations to emerge.
On the other hand, there is no doubt that genetics is a key player influencing concurrent changes between regions in morphometric, functional, and tissue properties (Alexander-Bloch et al., 2013;Brown et al., 2011;Goryawala et al., 2015;Shu et al., 2015;Yao et al., 2015).
Twin studies have demonstrated that specific anatomical co-variations rest on shared genetic influences (Alexander-Bloch et al., Schmitt et al., 2008Schmitt et al., , 2009Schmitt et al., , 2010. Myelination processes could be modulated by shared genetic influences due to a single gene, as a putative mechanism of synchronization of these processes among brain regions during neurodevelopment and aging. This mechanism has previously been proposed studying other brain anatomy variables (Meyer-Lindenberg et al., 2007;Pezawas et al., 2005). Along the same line of thought, the interaction of different genes involved in brain myelination (related to myelin structure, composition, development, or maintenance; Jungerius et al., 2008) may contribute to correlations in myelination between regions.
Further experimental work (animals and humans, in vivo and ex vivo) and theoretical modelling need to be conducted to support or disprove the impact of these putative factors and mechanisms on correlations in myelination between gray matter regions.

| Concurrent myelin changes in the brain: Aging effects
The study of myelination processes is crucial to understand brain anatomy and function. Here we extended the standard analysis of exploring myelin distributions across the gray matter (first-order analysis) to a second order approach, where the amount of myelin in each anatomical region is no longer relevant (the mean regional myelination values were detrended out before analysis). It is this approach that allows the estimation of the interactions between structures by the Pearson correlation coefficients. It has been proven that such type of analysis reveals hidden organizational properties of brain anatomy that are difficult to detect through classical univariate approaches (Bassett et al., light of the known increase in symmetry in brain function and anatomy with age expressed by the recruitment of homologous areas, which may operate to compensate for cognitive decline (Cabeza et al., 1997;Cabeza, 2002;Grady, McIntosh, Rajah, Beig, & Craik, 1999;Park and Reuter-Lorenz, 2009).
We also observed that the highest correlations were precisely between homologous regions for both age groups. This finding has been systematically reported in the literature using different functional and anatomical variables (Eyler et al., 2011;Mechelli et al., 2005;Melie-García et al., 2013;Sanabria-Diaz, Martínez-Montes, & Melie-Garcia, 2013;Schmitt et al., 2009). We hypothesize that genetic factors play an important role in fine-tuning the correlation of the myelination processes between homologous structures.

| Most connected structures in Myelin-Net, Myelin-Nets hubs, and aging effects
To better understand the topological organization of spatial myelination changes, we identified the structures where myelination processes are the most correlated with the myelination variations across the rest of the brain. Some structures were common to both age groups (e.g., bilateral supplementary motor cortex (SMC.L, SMC.R) and left superior parietal lobe (SPL.L)). The presence of SMC as a network hub has been described previously at different ages. In the case of the elderly connectome (76-94 years), this structure is preserved as a hub (Hwang, Hallquist, & Luna, 2013;Perry et al., 2015).
We evidenced a redistribution of hubs with aging ( Figure 7). There

| Organizational properties and aging effects of myelin nets
The topological properties of the Myelin-Nets were also modulated by aging. The "characteristic path length" in Old-Age was statistically higher than the Young-Age group. This alteration has been reported consistently in previous studies using other neuroimaging modalities (Gong et al., 2009b;Petti et al., 2016;Zhao et al., 2015;Zhu et al., 2012). The nature of this effect is unknown but could be related to compensatory effects and common pathological processes or shared vulnerability. These processes could influence the redistribution of principal myelin covariations (the highest correlation values) to specific regions, which ultimately may induce longer characteristic path lengths.
Additionally, we provide empirical evidences that after 48 years old this network attribute present dynamic changes until 75 years old.
The "clustering index" was increased in the Old-Age group. This is a measure of the similarity of myelin covariations among brain structure neighborhoods (not reduced to a physical neighborhood concept).
According to the graph theory, this increase could be generated by the establishment of new densely connected local clusters which may generate an uncontrolled "flow of information" through the entire network.
This measure is related to the local efficiency of the "information flow" of the networks and its abnormal performance could be attributed to shared vulnerability and compensatory mechanisms. The increased intralobe myelin covariation in 4 out of 7 lobes might explain the larger clustering index in old age.
The increasing "local efficiency" was accompanied by a "global efficiency" decline in Old-Age. In terms of the graph theory, these changes affect network performance pointing to a higher "wiring cost" for parallel "information transfer" between anatomical regions. The weakening of the global efficiency is explained by the larger characteristic path length present in old age.
The aging trajectory in both "clustering index" and "local efficiency" showed a minimum at 48 years old that is related to global myelin correlation strength changes and the cause producing that this attribute has a minimum value at the same age. The maxima of these two network properties and the maximum/minimum in characteristic-pathlength/global-efficiency are synchronized in some extend revealing an unbalance between segregation and integration processes in Myelin-Nets.
In our "Targeted Attack" analysis, we attempted to further investigate the age-associated changes in the resilience of the Myelin-Nets.
Networks were "attacked" at key nodes of the highest betweenness centrality (network hubs). We found that the Old-Age group showed higher resilience to these simulated attacks. A possible explanation for this finding may rest on the observed increase in symmetry and the presence of aberrant local circuits leading to increased "local efficiency" with age. Also the redistribution and the higher number of hubs in the Old-Age group favored a higher resilience to target damage. Attacks to hub regions in the Young-Age group produced devastating consequences in the integrity and stability of the network, where regional specialization is well defined.
The comparison of the networks of myelin reported in this study and Networks of anatomical covariance may bring new light on the interplay between tissue myelination and morphological changes in the brain. Our results indicate that Myelin-Nets and Networks of anatomical covariance have common topological features. For both brain networks, correlations between homologous regions are the highest and are modulated by age. Also, the global network attributes-characteristic path length, clustering index, and local efficiency-share the same age modulation across networks (Chen et al., 2011;Wu et al., 2012;Zhu et al., 2012). Additionally, Myelin-Nets and Networks of anatomical covariance share common hub regions such as precuneus, middle and superior frontal gyrus, prefrontal cortex, posterior cingulate cortex, inferior temporal gyrus, and medial temporal gyrus (He et al., 2007;Wu et al., 2012). These similarities suggest possible common principles driving the topological organization of the synchronized myelination and morphological changes across the gray matter. However, a number of differences are also apparent which might be due in part to differences in the data and analysis methods used, which motivates further investigation to be conducted using the same experimental conditions (same subjects sample, brain parcellation, weighted/binary graphs, etc.) to ease the comparability of the results.
Finally, a comparison between Myelin-Nets and white matter structural networks (WM-Nets) would shed light on the common principles driving the correlated myelination processes and anatomical (axonal) connectivity. A separate study, preferably based on Myelin-and WM-Nets extracted from the same cohort, is necessary to highlight the interplay between these processes. However, some similarities and differences can be highlighted in light of the main anatomical connectivity findings reported in the literature. First, in both myelin and anatomical networks, the precuneus and posterior cingulate gyrus are observed as centrally connected regions, independently of age (Gong et al., 2009a, b;Iturria-Medina, Sotero, Canales-Rodríguez, Alem an-G omez, & Melie-García, 2008). This observation is consistent with the study of (Hagmann et al., 2008), who identified a structural core within posterior medial and parietal cortex in the cortical anatomical network. A region equivalent to the precuneus was also observed as a hub in the macaque cortical network (Sporns, Honey, & K€ otter, 2007). Also, putamen and superior parietal structures were identified as most vulnerable areas in WM-Nets (Iturria-Medina et al., 2008) as we found in Myelin-Nets.
The topological efficiency of the WM-Nets exhibits an inverted Ushaped trajectory across the lifespan, peaking around the third decade of age (Zhao et al., 2015). Similarly, we found that the local efficiency of the Myelin-Nets followed an inverted U-shaped between the ages of 48 and 75 years old, but with a peak between 60 and 65 years old. In contrast, the global efficiency of the Myelin-Nets was found to follow a U-shape with a minimum between 55 and 60 years. While the clustering index and characteristic path length of WM-Nets were found to follow an inverted U-shaped trajectory (Zhao et al., 2015), these topological attributes were found to exhibit different aging trajectories in the Myelin-Nets in the age range 48-75 years. Interestingly, while the WM-Nets become less connected with age (Gong et al., 2009b), the connectivity in Myelin-Nets shows a tendency to increase with age from 48 years old onward.

| Further considerations and future work
Some issues should be addressed in future works. We propose the use of parcellations closer to the myelo-architectonic organization of the cortex to study Myelin-Nets topological properties. The Neuromorphometrics, AAL, Brainnetome, and Gordon atlases employed in this study do not take into account this key organization feature of the cerebral cortex. Partial correlation, instead of the Pearson correlation, should be favored to compute the interaction between anatomical structures accounting for the effect of other structures and possible global variables. An increase in the number of subjects per age group in an independent dataset would provide evidence for the reproducibility of our results.
A natural extension of the presented work is the parallel characterization of morphological and myelin covariance networks in gray matter. This combined study, combining micro and mesoscopic measures of the brain, would shed a new light into the biophysical mechanisms underlying the emergence of such networks. It should be highlighted that myelin is the primary contributor to MR image contrast (Geyer, Weiss, Reimann, Lohmann, & Turner, 2011), with a clear impact on local gray matter volume estimates Lorio et al., 2016). The parallel study of morphological and myelination change would therefore allow disambiguating true morphological brain changes and spurious changes arising from changes in myelination.
An important line of future research is the exploration of potential deviations in the correlation of myelination processes in pathologies like Alzheimer, Schizophrenia, Multiple Sclerosis, and Epilepsy. Furthermore, the combined study of anatomical, myelin, functional brain networks, by combining structural DWI, magnetization transfer-MRI, and MELIE-GARCIA ET AL. | 1549 functional MRI techniques in a same set of subjects would be tremendously useful for discovering similarities and differences to brain network properties obtained from different physiological variables. This would help to understand how changes in the topological organization of myelination processes are related to brain function. The study of the association between gray and white matter myelination processes over the lifespan will help to improve our knowledge about the basic principles of the brain topological organization. Finally, we find two potential limitations of our methodology: (a) some of the identified age-related differences in Myelin-Net topological organization may be because of possible registration errors in the Neuromorphometrics atlasing labeling methodology. However, these sources of bias were minimized by using a nonlinear advanced registration algorithm (geodesic shooting registration; Ashburner and Friston, 2011) and a probabilistic labeling approach. The influence of our results by the utilization of surfacebased and volume-based registration algorithms should be evaluated in future studies. (b) Although myelination largely dominates the MT contrast other processes such as inflammation (B elanger, Allaman, & Magistretti, 2011;Gloor, Scheffler, & Bieri, 2008), metabolism (Giulietti et al., 2012), and pH changes (Gillies, Raghunand, Garcia-Martin, & Gatenby, 2004;Gloor et al., 2008;Kucharczyk, Macdonald, Stanisz, & Henkelman, 1994) appear to contribute to this signal (Harrison et al., 2015). The neuroinflammation, for instance, is present during aging (Chung et al., 2009;Franceschi et al., 2007); therefore, our results may be influenced by the heterogeneity of this process in our age groups.

| C ONC LUSI ON S
In this article, we provide empirical evidence that myelination processes are spatially correlated across cortical gray matter. The patterns of myelination covariance show specific topological organization revealed using a graph theoretical approach. A number of specific correlated myelination phenomena were strongly modulated by age. In particular, we brought evidence that aging increases synchronicity between homologous regions. Our results are an important step toward elucidating the organizational principles behind the dynamics of the human brain anatomy across the lifespan. In particular, myelination processes are at the cross-way of several psychiatric and neurodegenerative brain diseases and are therefore crucial for their understanding.