|Home | About | Journals | Submit | Contact Us | Français|
Research in the macaque monkey suggests that cortical areas with similar microstructure are more likely to be connected. Here, we examine this link in the human cerebral cortex using 2 magnetic resonance imaging (MRI) measures: quantitative T1 maps, which are sensitive to intracortical myelin content and provide an in vivo proxy for cortical microstructure, and resting-state functional connectivity. Using ultrahigh-resolution MRI at 7 T and dedicated image processing tools, we demonstrate a systematic relationship between T1-based intracortical myelin content and functional connectivity. This effect is independent of the proximity of areas. We employ nonlinear dimensionality reduction to characterize connectivity components and identify specific aspects of functional connectivity that are linked to myelin content. Our results reveal a consistent spatial pattern throughout different analytic approaches. While functional connectivity and myelin content are closely linked in unimodal areas, the correspondence is lower in transmodal areas, especially in posteromedial cortex and the angular gyrus. Our findings are in agreement with comprehensive reports linking histologically assessed microstructure and connectivity in different mammalian species and extend them to the human cerebral cortex in vivo.
Comprehensive research in macaque monkeys has revealed an intricate link between cortical microstructure and connectivity (see Pandya et al. 2015 for a review). In particular, long-range connections preferentially occur between areas with similar microstructural features (e.g. Pandya and Sanides 1973; Pandya and Yeterian 1985, 1990; Barbas and Pandya 1989). Microstructural similarity has therefore been suggested as a general connectivity principle—or wiring rule—of the mammalian cortex (Barbas 2015; Pandya et al. 2015). This claim is supported by quantitative analyses of the macaque (Beul et al. 2015), cat (Beul et al. 2014), and mouse (Goulas et al. 2016a) connectomes, but remains largely unexplored in the human brain.
The aforementioned animal studies rely on combined tract-tracing and postmortem histology, invasive approaches that are unsuitable to study the human brain. While magnetic resonance imaging (MRI)-based techniques to assess cortical connectivity have been widely adopted, the analysis of human cortical microstructure remains challenging. Some recent studies have related cytoarchitectonic data, available from the literature, to connectivity measures derived from MRI in a separate group of subjects. Thus, van den Heuvel et al. (2015, 2016) demonstrated that the degree of connectivity of an area is related to the size and complexity of its layer III pyramidal neurons, a finding that is in line with previous work in macaque monkeys (Scholtens et al. 2014). Moreover, Goulas et al. (2016b) showed that the existence of connections between 2 areas is partly predicted by their similarity in supragranular cell density. These findings lend strong support for a relationship between macroscale connectivity and microstructural features in the human cortex. Here, we investigate this link in vivo, using high-resolution MRI measures, sensitive to connectivity and microstructural features, both acquired from the same subjects.
With recent advances in MRI technology, it has become possible to derive measures related to cortical microstructure in the living human brain. Classic histological approaches to describe cortical microstructure focus on the distribution of cells (“cytoarchitectonics” e.g. Brodmann 1909) or myelinated fibers (“myeloarchitectonics” e.g. Vogt and Vogt 1919) within the cortical sheet. Because the longitudinal relaxation time (T1) in MRI is sensitive to gray-matter myelin content (Bock et al. 2009; Stüber et al. 2014), maps of intracortical T1 have been introduced as an in vivo proxy for cortical microstructure and revived interest in myeloarchitectonic approaches (Nieuwenhuys 2013). These maps exhibit a decrease in estimated myelin content from primary to transmodal areas (Glasser and Van Essen 2011; Tardif et al. 2015), which is in line with findings from histological studies (e.g., Hopf 1955, 1956; Hopf and Vitzthum 1957). The co-occurrence of local changes in intracortical T1 maps with certain architectonic (Geyer et al. 2011; Glasser and Van Essen 2011), functional (Bridge et al. 2005; Sigalovsky et al. 2006), and topographic (Dick et al. 2012; Sereno et al. 2013) boundaries, as well as rapid changes in functional connectivity patterns (Glasser et al. 2016), exemplifies their usefulness to investigate cortical organization.
Cortical connectivity patterns are now commonly investigated in vivo through resting-state functional connectivity (Biswal et al. 1995). While functional connectivity can reflect indirect links between areas (Adachi et al. 2012), it is largely constrained by anatomical connections (Skudlarski et al. 2008; Honey et al. 2009; Hermundstad et al. 2013).
With this study, we investigate the relationship between in vivo measures sensitive to intracortical myelin and functional connectivity in the human cerebral cortex. We employ ultrahigh field MRI and dedicated image processing tools to compare functional connectivity patterns with cortical differences in T1-based myelin estimates at a submillimeter resolution. Motivated by the coherent findings in other mammalian species, we hypothesize that similarity in intracortical myelin relates to functional connectivity, thereby constituting a wiring rule of human cerebral cortex.
MRI data were acquired and published previously in 2 independent studies (Gorgolewski et al. 2015; Tardif et al. 2016). Here, we use data sets of 9 subjects who took part in both studies (5 females, mean age = 24.8 years, standard deviation = 1.2 years, see Supplementary Table 1 for subject identifiers). All subjects gave written informed consent and the studies were approved by the Ethics Committee of the University of Leipzig. Data were acquired on a 7T MR scanner (MAGNETOM, Siemens Healthcare) equipped with a 24-channel Nova head coil. In the first study (Tardif et al. 2016), a structural scan was acquired at an isotropic resolution of 0.5 mm using the MP2RAGE sequence (Marques et al. 2010, inversion time (TI)1/TI2 = 900/2750 ms, repetition time (TR) = 5 s, echo time (TE) = 2.45 ms, α1/α2 = 5°/3°, bandwidth = 250 Hz/px, echo spacing = 6.8 ms, partial Fourier = 6/8, scan time = 28 min per hemisphere) and an optimized radiofrequency pulse to minimize sensitivity to B1 inhomogeneity (Hurley et al. 2010). The resulting images were used to generate whole brain T1 maps. In the second study (Gorgolewski et al. 2015), a total of 4 whole brain resting-state scans were acquired in 2 sessions at an isotropic resolution of 1.5 mm, using a T2*-weighted echo-planar imaging sequence (field of view = 192 × 192 mm2, 70 slices, TR = 3 s, TE = 17 ms, α = 70°, bandwidth = 1116 Hz/Px, partial Fourier 6/8, GRAPPA acceleration with iPAT factor of 3, 300 volumes per scan). In the first of these 2 sessions, an additional short structural scan (11 min) was acquired using the MP2RAGE sequence with similar parameters as above but an isotropic resolution of 0.7 mm. The uniform T1-weighted (T1w) image generated from this scan was used as an intermediate registration target to project the resting-state data into the space of the high-resolution T1 map and will be referred to as the “low-resolution T1w image.”
All data processing and analyses were performed on Linux servers running Ubuntu 12.04.5 LTS. The following software tools were used for data processing: CBS Tools (v3.0, Bazin et al. 2014) as a plugin for MIPAV (v7.0.1, McAuliffe et al. 2001) and JIST (v2.0, Lucas et al. 2010), ANTs (v2.1.0, Avants et al. 2011), Nipype (v0.11.0, Gorgolewski et al. 2011), Nilearn (v0.2.3, Abraham et al. 2014), Meshlab (v1.3.0, meshlab.sourceforge.net), Freesurfer (v5.3.0, Dale et al. 1999; Fischl et al. 1999), and AFNI (v16.1.28, Cox 1996). Custom scripts were written in Python 2.7 using numerous libraries (Jones et al. 2001; McKinney 2010; Seabold and Perktold 2010; Pedregosa et al. 2011; van der Walt et al. 2011; Droettboom et al. 2016; Waskom et al. 2016). More details and all code can be found at https://github.com/juhuntenburg/myelinconnect.
Preprocessing of the high-resolution T1 maps was carried out using CBS Tools (Fig. (Fig.1).1). Tissue segmentation and reconstruction of white matter (WM), midcortical, and pial surfaces was performed as described previously (Bazin et al. 2014). Cortical thickness between the WM and pial surface was calculated. For each subject, 9 intracortical surfaces were generated between the WM and pial surface. The radial position of each surface was determined using a volume-preserving approach (Waehnert et al. 2014), which accounts for the dependence of layer thickness on cortical folding described by Bok (1929). For each node of the subject's midcortical surface, T1 values were sampled in radial direction on the surfaces (WM, 9 intracortical, pial), yielding a vector of 11 values that was assigned to each node. We will refer to this vector as the “T1 profile” (cf. Waehnert et al. 2014, 2016; Dinse et al. 2015). We used the multimodal, multicontrast surface registration approach (MMSR), introduced by Tardif et al. (2015), to derive a study-specific surface template from the high-resolution T1 maps. The resulting group-average surface mesh was simplified from approximately 700k to 70k nodes using MeshLab's Quadric Edge Collapse Decimation filter (percent reduction = 0.11, quality threshold = 1.0, preserving mesh boundaries, normals, and topology, optimal placement, postsimplification clean). This template surface provides the common reference space for subsequent analyses. For visualization purposes, the group template surface was partly inflated using Meshlab's Taubin Smooth filter (lambda = 0.6, mu = −0.3, steps = 150). For coregistration of the resting-state data, the low-resolution T1w images were preprocessed using FreeSurfer's standard recon-all pipeline.
Preprocessing of the resting-state time series was streamlined in a reusable pipeline using Nipype. For each resting-state scan, the first 5 volumes were removed to ensure steady-state magnetization. After simultaneous slice time and motion correction, as implemented in Nipy (Roche 2011), a temporal median image was created and bias field corrected using N4 as implemented in ANTs (Tustison et al. 2010). This median image was coregistered to the high-resolution T1 map of the same individual in 2 steps: in a linear step, both the median image and the high-resolution T1 map were registered to the low-resolution T1w image using Freesurfer's boundary-based registration with 6 degrees of freedom (Greve and Fischl 2009). The 2 resulting transformations were concatenated and used to project the median image into the space of the high-resolution T1 map. After carefully masking both images, a second nonlinear registration step was performed using ANTs’ symmetric diffeomorphic transformation model (SyN) and fast cross-correlation similarity metric (Avants et al. 2008). This last step is intended to account for nonlinear distortions in the functional data. All transformations were combined and used to project data from the individual's functional to the high-resolution structural space in a single interpolation (Fig. (Fig.1).1). Voxel-wise temporal denoising of the slice time and motion corrected time series was performed using Nilearn. Confounds that were removed from the time series included linear trends, the 6 motion regressors and their first derivatives, intensity outliers (z > 3) and motion outliers (mean composite norm > 0.5 mm), as identified by Nipype's rapidart algorithm, as well as regressors reflecting the signal in the WM and cerebrospinal fluid (CSF), following the aCompCor approach (Behzadi et al. 2007). For the latter, WM and CSF masks were derived using the MGDM segmentation as implemented in CBS Tools (Bazin et al. 2014), and eroded by 4 (WM) or 2 (CSF) voxels, respectively, to avoid partial volume effects with the gray matter. Nilearn's high_variance_confounds method with default settings was used to extract 5 high variance temporal components from voxels within each of the 2 masks, resulting in a total of 10 components that were included in the regression. The denoised time series were temporally normalized and bandpass filtered to a frequency range of 0.01–0.1 Hz. One subject was excluded due to excessive motion (composite norm max > 3 mm, mean > 0.15 mm across the 4 resting-state scans) and insufficient coregistration quality as assessed by visual inspection (Supplementary Fig. 1). Data of the remaining 8 subjects passed these criteria and were used for subsequent analyses.
The medial wall was delineated by hand on both hemispheres of the inflated group template surface mesh using AFNI's Surface Mapper. Additionally, masks were created to exclude nodes with signal quality of at least 1.5 standard deviations below the mean in at least one session of one subject and one imaging modality. Signal quality in each subject's structural image was derived from the second inversion images of the MP2RAGE sequence: the probability that a data point represents foreground (data, modeled by a uniform distribution) rather than background (noise, modeled by an exponential distribution) was calculated based on a distribution mixture model. The derived measure represents a heuristic definition of available signal in the T1 maps. Signal quality in the resting-state images was assessed separately for each subject and session using the temporal signal-to-noise ratio (tSNR). After slice time and motion correction, but before nuissance regression, each voxel's tSNR was calculated by dividing its temporal mean BOLD signal by the temporal standard deviation of its BOLD signal. All masks were combined and used to exclude low fidelity regions from subsequent analyses.
For each subject, the 5 central values (45%) of the T1 profiles were averaged to create a single intracortical average T1 value at each surface node, which is minimally biased by partial volume effects with the WM and CSF (Figs (Figs11 and and2).2). Throughout the manuscript, we assume this intracortical average T1 value to largely reflect intracortical myelin content. Limitations of this assumption will be addressed in “Discussion” section. Intracortical average T1 values were slightly smoothed on the individual midcortical surface using a Gaussian kernel with full-width-half-maximum (FWHM) of 1.5 mm, to prepare them for sampling on the lower resolution template surface. The group template surface was projected into each individual's native space using the MMSR-derived transformations and the intracortical average T1 was sampled from the closest point on the subject's original high-resolution surface. Thus, each node on the group-level surface was associated with a single intracortical T1 value for each subject. These values were subsequently averaged across subjects to derive the group-level intracortical T1 map. A T1 difference matrix was generated from this map by calculating the absolute difference in group-level intracortical T1 for each pair of nodes on the surface template. High values in this matrix therefore indicate node pairs that differ substantially in their intracortical average T1 values while values close to zero indicate pairs with very similar values. (See Fig. Fig.22 for a schematic of the processing steps..
The fully preprocessed resting-state time series were projected into the individual high-resolution structural space. Subsequently, the time series data of each subject were sampled onto the group template surface at midcortical depth and smoothed along the cortex using a Gaussian kernel (FWHM = 3 mm). Due to its lower spatial resolution, the resting-state data were not sampled at different intracortical depth levels. Whole brain functional connectivity matrices were derived for each subject and session by calculating the Pearson's product–moment correlation between the time series of each pair of nodes. All matrices were Fisher r-to-z-transformed and averaged across subjects and sessions. The resulting average functional connectivity matrix was back transformed to Pearson's r values. (See Fig. Fig.22 for a schematic of the processing steps.)
Pearson's product–moment correlation coefficient was calculated between the upper triangles of the group-average functional connectivity and T1 difference matrices. Moreover, functional connectivity and T1 difference patterns were related to each other in a node-wise fashion: Pearson's product–moment correlation coefficient was calculated between each row of the functional connectivity matrix (representing functional connectivity of one node to all other cortical nodes) with the same row in the T1 difference matrix (representing the T1 difference of the same node to all other cortical nodes). This results in one value per node, indicating how well its functional connectivity to any other node can be predicted from its T1 difference to that node (and vice versa). The resulting map was compared with the group average of signal quality measures in both imaging modalities, assessed as described in “Masks” section (see Supplementary Fig. 2). To control for the effect of spatial distance, the above analyses were repeated after regressing Euclidean distances between nodes in 3D space against both functional connectivity and T1 difference. Thus, both the upper triangle correlation and node-wise correlations were recalculated for the residuals of the functional connectivity and T1 difference after distance regression. Euclidean distance was also correlated to T1 difference and functional connectivity separately (Supplementary Fig. 3). Importantly, Euclidean distance is a biased estimate of the actual length of fiber tracts through the WM, which systematically underestimates certain pathways. Inferences drawn using it as a measure of connection length must therefore be evaluated in the light of this serious limitation. To control for the relationship between cortical thickness (CT) and intracortical T1, CT was smoothed (FWHM = 1.5 mm), sampled in the subjects’ native space and averaged on the group template surface. A matrix of absolute difference in group-level CT between each pair of nodes was regressed against the T1 difference matrix, and vice versa, before correlation to functional connectivity.
To disentangle which particular functional connectivity features underlie the relationship to intracortical T1, we decomposed the functional connectivity matrix into a set of one-dimensional components, reflecting distinct aspects of cortical connectivity (see Fig. Fig.2).2). For this, the group average functional connectivity matrix was transformed into a positive similarity matrix L through increasing each element by 1 and then dividing it by 2. Thus, nodes with highly correlated signals are considered similar, while nodes with highly anticorrelated signals are considered dissimilar. Nonlinear dimensionality reduction of L was performed using diffusion maps (Coifman and Lafon 2006, Implemented in https://github.com/satra/mapalign.) In short, diffusion maps embed the high-dimensional connectivity data in low-dimensional space by first transforming the similarity matrix L in a Markov chain with transition probability matrix M:
If the diffusion operator α is set to 1, M is equal to the classical normalized graph Laplacian. Here, we use α = 0.5 in which case transition in the Markov chain approximates Fokker–Planck diffusion and the method is less sensitive to nonuniform sampling of the data on the underlying manifold. The eigenvectors of M constitute the coordinate system of the embedded space. Here, we embed the 70k × 70k similarity matrix L in 100 dimensions, which we will refer to as functional connectivity components FC1–100. These components are naturally ordered by decreasing eigenvalues, that is, FC1 explains most of the variance, etc. To approximate the amount of variance explained by individual components, we make the simplifying assumption that FC1–100 together explain 100% of the variance (FC100 explains 0.29% of this variance). Note that these values refer to variance in the transition matrix M, not in the original functional connectivity matrix.
Diffusion maps have been used for the analysis of resting-state fMRI data before and in-depth discussions can be found in Langs et al. (2015a, 2015b) and Margulies et al. (2016). In contrast to linear decomposition approaches, such as principal component analysis, nonlinear methods retain both local and global relationships between data points in the embedded space without requiring kernel manipulations. Diffusion maps, specifically, are more robust to noise in the connectivity matrix than other nonlinear dimensionality reduction techniques, such as Isomap (Tenenbaum et al. 2000).
Pearson's product–moment correlation coefficient was calculated between the first component of the functional connectivity embedding (FC1) and intracortical T1. To account for heteroscedasticity, we assessed the significance of this correlation using the robust sandwich variance estimator (White 1980) as implemented in Statsmodels. Higher order functions were fitted as well, to assess whether they would provide a better description of the relationship between FC1 and T1. In a next step, intracortical T1 was modeled as a linear combination of one or multiple FC components using Ordinary Least Squares regression as implemented in Scikit-learn such that:
We compared all possible combinations of the first 20 FC components (FC1–20, each explaining >1%, together explaining 61.6% of the variance in M, 1 048 575 different models). Again, significance of the coefficients was tested using the robust sandwich variance estimator.
To compare the performance of combinations of FC components to model intracortical T1, we derived the Bayesian information criterion (BIC) as:
where p is the number of parameters and n the number of data points used for fitting, R is the data range, is the residual variance and the sum of squared residuals. The derivation of the BIC closely follows that in Bazin and Vezien (2005) and can be found in the Supplementary Methods. We estimated from local variance in the group average T1 data using all surface nodes that were not excluded by the mask:
where is an estimate of the median absolute difference between neighboring T1 values. Assuming that these differences follow a half Gaussian distribution their median is given by and . The data range R is set to the range of group average T1 values not excluded by the mask.
Random data sets were simulated through the following steps: first, values for each node on the surface template were drawn from a normal distribution, creating 1000 data sets per hemisphere. Second, data were smoothed on the surface applying differently sized Gaussian smoothing kernels (FWHM = 1.5, 3, 6, 12, and 24 mm). Third, after masking, the values of each random data set were adjusted to those of the original intracortical T1 map using histogram matching. Each random data set was fitted as described for the T1 data in Functional Connectivity Decomposition and Comparing Components to T1, employing the models consisting of either FC1 alone or a linear combination of FC1, 5, and 6 (the best performing model in the previous analysis). The coefficient of determination was calculated for all models and used as null distribution to assess significance of the fit achieved for the real intracortical T1 map. Note that for this analysis, the left and right hemispheres were treated separately for both the random and the real T1 data. Using whole brain data would result in an unfair comparison because both the FC components and the original T1 map show high interhemispheric symmetry while the random data do not.
We created 2 group average matrices: one containing the difference in intracortical T1 between each pair of nodes on the surface (“Intracortical T1 and T1 difference matrix”) and a second one with the functional connectivity between each pair of nodes (“Functional connectivity matrix and comparison to T1 difference matrix,” see also Fig. Fig.2).2). We found both matrices to be anticorrelated (Pearson's r = −0.34, P < 0.0001, ). Thus, the more similar 2 nodes are in their T1 values, the higher their functional connectivity.
To investigate whether this inverse relationship between functional connectivity and T1 difference is homogeneous across the cortex, we also calculated the correlation in a node-wise fashion. For each surface node with sufficient signal quality in both imaging modalities (cf. “Masks”), we compared its functional connectivity with all other nodes with the respective T1 differences to those nodes. This results in one value per surface node, indicating how well T1 difference predicts functional connectivity for this particular node. Figure Figure33 shows that the relationship between functional connectivity and T1 difference follows a distinct pattern across the cortex. Unimodal areas display a particularly strong link between the 2 measures. Conversely, regions that show little correlation between functional connectivity and T1 differences include posteromedial cortex, superior portions of the inferior parietal lobule, superior temporal sulcus and middle temporal gyrus, medial prefrontal regions, anterior insular cortex and posterior portions of the superior and middle frontal gyrus. The correspondence between intracortical T1 and connectivity strength is thus generally stronger in unimodal than in transmodal areas. This effect shows little dependence on the signal quality in the structural (Pearson's r = −0.18, P < 0.0001) and resting-state (Pearson's r = −0.09, P < 0.0001) data. In fact, particularly low correspondence of functional connectivity and T1 difference in posteromedial cortex and the inferior parietal lobule coincides with high signal quality in both imaging modalities (Supplementary Fig. 2).
We observed the same general pattern of correspondence between functional connectivity and differences in a T1w/T2w-based estimate of intracortical myelin content in a large, independent data set (n = 820), made available through the Human Connectome Project (HCP, Glasser et al. 2013; van Essen et al. 2013, see Supplementary Methods and Supplementary Fig. 7 for details).
It is well known that connectivity decreases with spatial distance (e.g., Young 1992; Scannell et al. 1995) and reasonable to assume a similar relationship for microstructural similarity. We therefore investigated whether the observed relationship between functional connectivity and T1 difference is driven by a common dependence on spatial distance. To this end we created a third matrix, describing the Euclidean distance of each pair of nodes on the surface template. Of note, Euclidean distance is a biased estimate of connection length and the ensuing results must be interpreted in due consideration of this limitation. As expected, functional connectivity and Euclidean distance were highly correlated (Pearson's r = −0.41, P < 0.0001, , Supplementary Fig. 3a ). However, T1 difference was independent of Euclidean distance when calculated across the whole cortex (Pearson's r = 0.01, P < 0.0001, Supplementary Fig. 3b ). After regressing Euclidean distance against both T1 difference and functional connectivity, the anticorrelation between residual T1 difference and residual functional connectivity remained essentially unchanged (Pearson's r = −0.37, P < 0.0001, ).
We also repeated the node-wise analysis after regressing T1 difference and functional connectivity against Euclidean distance for each node separately. The overall pattern of the relationship between functional connectivity and T1 is not affected by distance regression (Supplementary Fig. 4). Only some areas, such as supplementary motor area, midcingulate cortex, and medial temporal gyrus, show increased anticorrelation. Thus, the relationship between functional connectivity and T1 is largely independent of a shared dependence on spatial distance.
Beul et al. (2015) found cortical thickness (CT) to be related to microstructural differentiation and cortical connectivity. However, in their study, the association between CT and connectivity disappeared when controlling for microstructural similarity, while microstructure was still related to connectivity when controlling for thickness. We found a similar pattern in our data. Difference in CT between nodes showed a strong correlation to T1 difference (Pearson's r = 0.70, P < 0.0001) and a moderate correlation to functional connectivity (Pearson's r = −0.28, P < 0.0001). Regressing CT difference against T1 difference reduced, but did not remove, the correlation between T1 difference and functional connectivity (Pearson's r = −0.20, P < 0.0001 vs. r= −0.34, P < 0.0001 before CT regression). However, CT difference was no longer related to functional connectivity after T1 difference had been regressed from it (Pearson's r = −0.06 , P < 0.0001). Intracortical T1 thus seems to capture variance in functional connectivity patterns beyond what is accounted for by CT alone.
After establishing an association between functional connectivity and T1 difference in high-dimensional space, we sought to identify the low-dimensional features of functional connectivity that specifically underlie this link. To this end, we decomposed the group-level functional connectivity matrix using nonlinear dimensionality reduction. The result are 100 components (FC1–100), each representing a distinct aspect of functional connectivity in a single dimension. The value assigned to a given node in each of these dimensions indicates its position along a spectrum of connectivity patterns: nodes with similar values on a given component resemble each other in their connectivity pattern, or more precisely, in the aspect of connectivity that is captured by this particular component. Nodes with very different values on the same component have more distinct connectivity profiles. Unlike independent component analysis, the method applied here does not enforce spatial independence of different components, but results in a set of gradients spanning the entire cortex (see “Functional Connectivity Decomposition and Comparing Components to T1” section for more details). This approach enables us to describe the relationship between the distribution of intracortical T1 and the main modes of variation of connectivity patterns across the cortex.
We initially focused on the principal gradient (FC1) that captures the main variation of functional connectivity patterns across the cortex (12% of the variance). As visible in Figure Figure4a,4a, FC1 spans a gradient between unimodal visual, motor, somatosensory, and auditory areas at one end of the spectrum, and higher order association areas in frontal, parietal, and temporal cortex at the other (see also Supplementary Fig. 5).
The one-dimensional representation of functional connectivity in FC1 can directly be compared with the group-average intracortical T1 map (see Fig. Fig.2).2). This map is used as an estimate of intracortical myelin content, with higher T1 values indicating less myelin (Waehnert et al. 2016). Limitations of this assumption will be discussed below. It reveals the highest estimated intracortical myelin content in primary visual, somatosensory, auditory, and motor regions. Intermediate values are found in posteromedial cortex, especially ventral portions adjacent to the corpus callosum, as well as superior parietal cortex, and relatively low myelin content prevails in prefrontal cortex, middle and inferior temporal gyrus as well as the inferior parietal lobule (Fig. (Fig.4a4a and Supplementary Fig. 5).
Node-wise comparison of the T1 and FC1 map reveals a substantial spatial correlation (Fig. (Fig.44b, Pearson's r = 0.53, P < 0.0001, ), exceeding the overall correlation in high-dimensional space. This implies that FC1 specifically captures an aspect of functional connectivity, which is related to intracortical T1. However, clear differences between the spatial layout of FC1 and T1 remain, most prominently in the superior part of the inferior parietal lobule and posteromedial cortex. Moreover, the bivariate distribution indicates a nonlinear relationship between FC1 and T1 (Fig. (Fig.44b). We therefore fitted higher order models to the data and compared them using an adapted version of the BIC (see “Model Comparison” section). Although more complex functions, such as polynomials or a sigmoid, explained the data better, their overall performance was worse, due to the higher number of parameters in these models (Supplementary Fig. 6).
Again, we found a similar relationship between FC1 and T1w/T2w-based intracortical myelin content in the HCP data set (n = 820, Pearson's r = 0.52 , P < 0.0001, see Supplementary Methods and Supplementary Fig. 8 for details).
The first component of the functional connectivity decomposition captures the main variance in the data. Still, a considerable amount of variance is explained by subsequent components (Fig. (Fig.55a), which might be essential to establish the link to intracortical T1. We therefore investigated whether a linear combination of multiple FC components can improve the fit to intracortical T1. We focused on components that each explains more than 1% of the variance (FC1–20). Every possible linear combination of these 20 components was fitted to the intracortical average T1 data. We then compared all models using the same implementation of BIC as in the previous section. A lower BIC value indicates a better balance of model fit and complexity.
Figure Figure55b shows the best performing model for each group of models employing the same number of FC components. Importantly, the model selection procedure is stable in that the best model in each group is equal to the best model in the previous group, except for the addition of one further. Moreover, the results are resistant against modest variation in the empirically determined noise prior (Supplementary Table 2). The most parsimonious model to fit intracortical T1 is a linear combination of component FC1, 5, and 6:
FC5 (explaining ca. 4% variance) contrasts posteromedial cortex and the superior parietal lobule at one end of the spectrum with inferior frontal gyrus, anterior insula, superior temporal cortex, ventral parts of precentral and postcentral gyrus, inferior parietal lobule, and visual cortex at the other (Fig. (Fig.55c, top). The most remarkable feature of FC6, which explains approximately 3% of variance, is that it separates primary sensorimotor regions from secondary regions and unimodal association areas within a given modality (Fig. (Fig.55c, bottom).
Combining FC1, 5, and 6 most notably improves the fit to intracortical T1 in posteromedial cortex (Fig. (Fig.66a). Moreover, visual regions become more similar in value to somatomotor and auditory regions than for FC1 alone. As shown in Figure Figure66b, including FC5 and FC6 also increases node-wise correlation between T1 and modeled T1 (Pearson's r = 0.67, P < 0.0001, , a quadratic model was fitted as well but without substantial gain, Supplementary Fig. 6). Differences to intracortical T1, however, remain. In particular, compared with the original T1 map, regions around the angular gyrus still display relatively high values in the combined FC1, 5, 6 map. The fitted T1 values also do not extend as far into the low range as original T1 values, resulting in residual differences especially in primary sensorimotor regions.
In order to ensure that the spatial correlation values obtained in the previous section do not result from spatial autocorrelation alone, we performed permutation testing. We simulated 1000 random data sets per hemisphere and smoothed them with different kernel sizes, ranging from the kernel applied to the original T1 maps (FWHM = 1.5 mm) to very large kernels (FWHM = 24 mm). We next fitted general linear models based on FC1 alone or FC1, 5, 6 to the random maps, as was done for the real data. Figure Figure77 shows that the model fit obtained for the real data far exceeds what is expected from data smoothness alone, even for large smoothing kernels (P < 0.005 in all cases). This is true for both hemispheres and both models, although the more complex model naturally fits the random data better than the simple model.
We demonstrated a systematic relationship between functional connectivity and a T1-based estimate of intracortical myelin in the human cortex. Regions with similar intracortical T1 show higher functional connectivity than regions that differ in intracortical T1. We initially illustrated the link to T1 in the high-dimensional space of the full functional connectivity matrix. By embedding functional connectivity into a set of one-dimensional components, we then identified particular functional connectivity features that relate to intracortical T1. Across both approaches, posteromedial cortex and superior portions of the inferior parietal lobule stood out as regions of low correspondence between functional connectivity and intracortical T1. In the following, we detail how the neurobiological signature of these regions, as well as characteristics of the MRI-based measures, can help to understand this effect.
In this study, we used intracortical T1 to estimate myelin content. Validation against histological data has confirmed that T1 contrast reflects myelin in cortical gray matter (Bock et al. 2009; Geyer et al. 2011; Stüber et al. 2014), displaying a gradient of decreasing myelin density from primary toward transmodal regions (cf. Vogt and Vogt 1919; Hopf 1955, 1956; Hopf and Vitzthum 1957; Sanides 1962). Iron concentration has also been shown to contribute to intracortical T1 contrast (Stüber et al. 2014) and an alternative interpretation of the current results is a systematic relationship between functional connectivity and iron levels. However, intracortical iron is strongly colocalized to myelin (Fukunaga et al. 2010), with ferritin particles embedded in the myelin sheath functioning as a storage for oligodendrocytes, that requires iron for the production and repair of myelin (Connor and Menzies 1996; Todorich et al. 2009). Therefore, independent of the exact contributions of iron and myelin to the T1 contrast, it appears justified to interpret T1 as to largely reflect the distribution of intracortical myelin.
We acquired quantitative T1 maps using the MP2RAGE sequence (Marques et al. 2010), like many recent works (Dick et al. 2012; Sereno et al. 2013; Tardif et al. 2015; Waehnert et al. 2016). This approach has the advantage of minimizing sensitivity to B1 inhomogeneities, allowing for quantitative intersubject and intersite comparison, and disentangling T1 from contribution other factors such as proton density and T2*, present in standard T1-weighted images (Turner 2015; Weiskopf et al. 2015). Alternative quantitative MRI techniques, such as magnetization transfer imaging (Dousset et al. 1992) and myelin water imaging (Mackay et al. 1994), are more specific to myelin than T1, but provide lower spatial resolution. Moreover, quantitative T1 has recently been characterized by the highest intrasubject and intersubject reliability in a comparison of several approaches to map intracortical myelin (Haast et al. 2016). Some studies have used a ratio of T1w over T2w images to assess intracortical myelin (Glasser and Van Essen 2011; Glasser et al. 2014, 2016). The relationship of T1w/T2w contrast to myelin content has not yet been validated and does not allow for quantitative comparison across studies. However, the spatial distribution of the T1w/T2w ratio generally appears similar to quantitative T1 maps and we were able to confirm the main results of the current study using T1w/T2w-based estimates of intracortical myelin. As myelin imaging becomes more abundantly used, it will be crucial to carefully compare different in vivo approaches with each other as well as with histological data.
For the in vivo assessment of cortical connectivity, we used resting-state functional connectivity, which has become a widely used tool in human neuroimaging. Large-scale functional connectivity patterns robustly reproduce across subjects (Damoiseaux et al. 2006), test–retest sessions (Shehzad et al. 2009; Zuo and Xing 2014), sites, and protocols (Biswal et al. 2010). Most importantly, resting-state functional connectivity is largely constrained by anatomical connections, as demonstrated in multiple studies in the human brain (Skudlarski et al. 2008; Honey et al. 2009; Hermundstad et al. 2013) and in comparison with tract-tracing data in monkeys (Vincent et al. 2007; Miranda-Dominguez et al. 2014). Since functional connectivity is based on temporal correlations, however, it does not necessarily reflect monosynaptic connections (Honey et al. 2009; Adachi et al. 2012). It is thus important to keep in mind that resting-state functional connectivity measures aspects of cortical organization that are closely related to, but not directly measuring, structural connections.
Microstructural similarity has been suggested as a candidate principle underlying the organization of cortical connections (Barbas 2015; Pandya et al. 2015). Respective links between microstructure and connectivity were established in the cortex of different mammalian species, by combining invasive tract-tracing and postmortem histology (Beul et al. 2014, 2015; Scholtens et al. 2014; Pandya et al. 2015; Goulas et al. 2016a). Due to the invasive nature of this approach, it cannot directly be translated to the human brain. Two recent studies have addressed this challenge by comparing diffusion-weighted MRI data to cytoarchitectonic atlas information (von Economo and Koskinas 1925), demonstrating a relationship between connectivity and layer III neuron complexity (van den Heuvel et al. 2015, 2016), as well as supragranular cell density (Goulas et al. 2016b). Spatial accuracy, however, is inevitably limited when mapping atlas information in stereotactic space and the discrete parcellation scheme prohibits the analysis of gradual changes or variable boundaries. We pursued an alternative approach, directly comparing in vivo measures related to cortical connectivity and microstructure in the same group of individuals at high spatial resolution. While these data facilitate the extraction of intracortical T1 values with high spatial precision, it is important to note that the current resolution does not allow the investigation of cortical features at the microscale. Differences in intracortical average T1 are likely to reflect differences in the underlying microstructure, namely the density of myelinated fibers. However, multiple microscale configurations can result in the same average T1 value when measured at a relatively coarse resolution of 0.5 mm and averaged across different cortical depth. Despite these limitations, we found that intracortical average T1 is related to functional connectivity, indicating that functional connectivity is higher between areas with similar intracortical myelin levels. Our results thus extend the aforementioned line of research to the living human brain and support the recognition of microstructural similarity as a general wiring rule of cerebral cortex.
Clearly, we do not mean to imply that microstructure is the sole determinant of connectivity. A well-established wiring rule states that 2 areas are more likely to be connected if they are close to each other in space (Young 1992; Scannell et al. 1995). This distance rule, however, cannot account for the existence of long-range connections and it must be assumed that additional factors are at play. In the current study, we confirm that a relationship exists between Euclidean distance and functional connectivity (Supplementary Fig. 3a ), but ensure that the correlation between functional connectivity and T1 is not driven by spatial proximity (Supplementary Fig. 4). While Euclidean distance is a biased estimate of connection length, entailing significant limitations of inferences drawn on its basis, our findings can be interpreted as a provisional indication that cortical microstructure and spatial distance are independent factors shaping functional connectivity (cf. Goulas et al. 2016a).
Moreover, similarity in cortical thickness (CT) has previously been shown to be related to connectivity in the macaque monkey (Beul et al. 2015) and we find it to be associated with functional connectivity in our sample. Similarly to Beul et al. (2015), we could show that this relationship disappears when controlling for the shared variance of CT and intracortical T1. In contrast, intracortical T1 captures variance in functional connectivity beyond what is explained by CT. Our results suggest that CT can be viewed as a pragmatic surrogate marker for cortical microstructure, capturing some of the same aspects of cortical organization as intracortical T1 (cf. Wagstyl et al. 2015).
While a general link between cortical microstructure and functional connectivity is supported by our data, we also find important deviations. Regions in which functional connectivity and intracortical T1 are least correlated consistently lie in transmodal cortex, especially in posteromedial cortex and the inferior parietal lobule. This pattern persists through different analytic approaches (Figs (Figs33 and and4)4) and cannot be attributed to the proximity of areas or the distribution of noise in the data (Supplementary Figs 2 and 4).
One possible explanation for the observed deviations are methodological challenges of our MRI-based approach. When visually comparing our T1-based measures of intracortical myelin with histologically derived maps of overall myelin density from Adolf Hopf (recently reviewed in Nieuwenhuys and Broere 2016), the general spatial pattern appears similar. Noteworthy differences, however, exist in posteromedial cortex and the inferior parietal lobule. In Hopf's map, these regions have lower myelin densities than most frontal regions, while the opposite relationship holds for our data (Fig. (Fig.44a) as well as for previously published myelin maps (e.g., Glasser and Van Essen 2011; Tardif et al. 2015). Strikingly, the same regions show the weakest link between functional connectivity and intracortical T1 in our study. It is therefore, possible, that the lower correspondence in posteromedial and inferior parietal regions is an artefact of in vivo derived T1-based myelin measures. Moreover, our group-average functional connectivity estimates could be impaired by interindividual differences in functional connectivity patterns, which are relatively high in inferior parietal regions (Mueller et al. 2013), although this does not apply to posteromedial cortex.
In addition to methodological differences, it is important to keep in mind that previous studies were performed in nonhuman primates or lower mammals. Compared with those species, transmodal areas in particular are thought to have undergone massive expansion in the human lineage (Hill et al. 2010). It has been suggested that this rapid expansion, by way of increasing the distance to molecular patterning centers, relaxed developmental constraints and allowed for new properties to emerge in transmodal areas (Buckner and Krienen 2013). Indeed, transmodal, as compared with unimodal, areas show denser and more variable connectivity patterns (Pandya et al. 1988; Mueller et al. 2013). Cortical folds in these areas emerge later during development (Hill et al. 2010), are more variable across individuals (Zilles et al. 1997), and show a less straightforward relationship to architectonically defined boundaries (Fischl et al. 2008). Myelination in transmodal areas also occurs late during development (Flechsig 1920) and remains relatively light in the adult brain. According to Braitenberg (1962), a major function of myelination in cortical gray matter is to inhibit the formation of new connections. Low myelin content in transmodal regions might therefore entail increased synaptic plasticity and make these areas more suitable to support adaptive behavior and learning. Finally, gene expression profiles distinguish resting-state networks (Krienen et al. 2016) and hubs (Vértes et al. 2016) related to sensorimotor functions from those in transmodal association and paralimbic regions.
Taken together, transmodal cortex appears to differ from unimodal areas in displaying more complex and variable connectivity and morphological patterns, as well as molecular and cellular characteristics that point toward increased plasticity. It could therefore be hypothesized that wiring in these regions is governed by flexible mechanisms as well. More rigid wiring rules, such as spatial proximity or microstructural similarity, might consequently have less impact. However, not all of transmodal cortex shows equally low correspondence between T1 and functional connectivity in our analyses. The most pronounced deviations occur in posteromedial cortex and the inferior parietal lobule, in particular the angular gyrus. Interestingly, these regions share a very distinctive connectivity profile. They have been identified as central hubs of structural connectivity in the human cortex (Hagmann et al. 2008) and, unlike most other regions, exhibit both high distant and high local functional connectivity (Sepulcre et al. 2010). Below we discuss how these connectional features might influence the correspondence of functional connectivity to intracortical T1 in posteromedial and inferior parietal cortex.
As described in the previous section, one characteristic of those areas in which functional connectivity and intracortical T1 are least associated is that they have particularly dense and diverse connectivity patterns. For instance, posteromedial cortex has heterogeneous functional and structural connectivity to unimodal, transmodal, and paralimbic areas (Margulies et al. 2009). It is conceivable that only a portion of these connections are explained by microstructural similarity. Only certain features of the functional connectivity patterns would then relate to intracortical T1, but the link is obscured when considering all functional connections of an area at once. We addressed this possibility by decomposing the functional connectivity space into different components, each representing specific aspects of functional connectivity as one dimension. We could thus describe the relationship between the distribution of intracortical T1 and the main modes of variation of connectivity patterns across the cortex.
The first component of this decomposition (FC1), spanning a gradient between primary sensorimotor and transmodal association regions (Fig. (Fig.44a, cf. Margulies et al. 2016), showed a stronger relationship to intracortical T1 than the full functional connectivity matrix. Deviations in transmodal regions, especially in posteromedial cortex, were still present in this comparison but substantially alleviated when including 2 additional components of the functional connectivity decomposition (FC 5 and 6) in a simple linear model of intracortical T1 (Figs (Figs55 and and6).6). Importantly, the added components were not the ones that explained most of the remaining variance, indicating that significant aspects of functional connectivity show no strong relationship to the distribution of intracortical T1. When scrutinizing those particular features which improve the fit to intracortical T1, it appears that FC6 captures fine differences between primary and higher order unimodal sensorimotor areas that are not present in FC1. The most important property of FC5 in the current context could be that it separates postcentral transmodal areas, which showed the strongest deviation from intracortical T1 when using FC1 alone, from transmodal areas in medial and ventrolateral (but not dorsolateral) frontal cortex. It should be noted that the foregoing interpretation of FC components remains tentative, as it is based on only a subset of salient features in the complex spatial patterns of these maps.
We consider decomposing functional connectivity patterns into overlapping gradients, rather than spatially discrete parcels, a promising new approach to investigate cortical organization. Recent work indicates that functional connectivity gradients closely relate to the functional specialization of the cortex (Haak et al. 2016; Margulies et al. 2016). The arrangement of cortical areas in spatial gradients has already been noted by Vogt and Vogt (1919) and Brockhaus (1940) and is fundamental to the evolutionary accounts of cortical organization by Sanides (1962) and Pandya et al. (2015). This concept does not conflict with the existence of areal boundaries, but rather attempts to situate individual areas within an overarching layout. Gradient-based approaches to describe cortical organization might therefore become a fruitful complement to discrete parcellations (e.g Glasser et al. 2016).
In this study, we used publicly available high quality MRI data (Gorgolewski et al. 2015; Tardif et al. 2016) and dedicated tools to harness the high spatial resolution (Bazin et al. 2014; Waehnert et al. 2014; Tardif et al. 2015). Nevertheless, our analyses and conclusions are subject to several limitations, which might partly be overcome in future studies.
One caveat is that the high field strength exacerbates the problem of MRI susceptibility artefacts. In our data, both structural and functional images suffer from low signal-to-noise ratio in orbitofrontal and ventral temporal areas. The functional data additionally display strong distortions in these areas, which cannot be corrected completely by nonlinear coregistration to the anatomy. To prevent unreliable signal in respective areas from driving our results, we decided to exclude them from our analyses (ca. 1.5% of surface nodes, see “Masks” section) and all results and statements only pertain to the included regions.
Further, given the distribution of our data (Figs (Figs44 and and6),6), it might be asked whether linear regression is the appropriate analytic approach. We modeled the data using higher order functions, which, however, did not entail substantial improvements when accounting for model complexity (Supplementary Fig. 6). To address the presence of heteroscedasticity, we used a robust sandwich variance estimator for significance testing (see “Functional Connectivity Decomposition and Comparing Components to T1” section). A more difficult question is, whether the multimodal distribution of the data reflects a principled division of the cortex in different classes which should be treated as discrete clusters. In our opinion, the considerable amount of data points falling between the denser data clouds indicate the existence of gradual changes in cortical features, which would be neglected by clustering approaches (see also “Investigating Connectivity Patterns Through Spatial Gradients” section). We thus favored a regression approach here, but a more in-depth treatment of this question may be subject to future studies.
We pursue a group-average approach using a study-specific template derived with a state-of-the art multimodal surface-based registration algorithm (Tardif et al. 2015). While such approaches have been shown to outperform volumetric and single modality registration techniques (Robinson et al. 2014; Tardif et al. 2015), intersubject averaging inevitably introduces inaccuracies in areas with high interindividual differences. Prospective work might circumvent this issue by investigating the relationship of functional connectivity and intracortical myelin within individual subjects (cf. Glasser et al. 2016).
Our estimate of intracortical myelin content takes into account the dependence of layer thickness on cortical morphology and minimizes partial volume effects. However, we use an intracortical average value, while classic myeloarchitectonic techniques are largely based on the differential distribution of myelinated fibers across the different cortical laminae (e.g., Vogt and Vogt 1919). Radial differences in myelin content can be comparable in size to differences between cortical areas (Lutti et al. 2014) and quantitative analyses of myelin or cell density profiles, running perpendicular to the cortical surface, have been used for observer-independent definition of areal boundaries in postmortem data (Schleicher et al. 1999; Annese et al. 2004; Eickhoff et al. 2005). Initial approaches to leverage the radial distribution of intracortical T1 in MR images in a similar way show promise that this information can be used to distinguish cortical areas (Dinse et al. 2015, Marques et al. 2016; Waehnert et al. 2016). Methods to extract and analyze such T1 profiles are still in their infancy, but more elaborate and robust approaches are likely to emerge from this line of research soon. Similarly, first studies demonstrating layer-resolved fMRI contrast (Goense et al. 2012; Hubert et al. 2015; Guidi et al. 2016) raise the possibility of differentiating functional connectivity between input and output layers. While the current study uncovers a general link between functional connectivity and intracortical T1, we expect that increasing spatial resolution in MRI data (Gallichan et al. 2015; Stucht et al. 2015) and improved image processing techniques will make it possible to refine this observation to the scale of individual laminar compartments in the near future.
The current resurgence of interest in human cortical anatomy is closely intertwined with methodological advances, providing the tools to investigate it in vivo. Some well-studied phenomena, emerging from decades of research in experimental animals, can now be explored noninvasively in the human brain. In the current study, we focus on a link between cortical microstructure and connectivity, which has consistently been described in animal studies, but until now had not been investigated in the living human brain. We demonstrate a systematic relationship between T1-based intracortical myelin content and resting-state functional connectivity. Our results support the view that microstructural similarity represents a general wiring rule of mammalian cortex. Although methodological difficulties persist and questions about underlying mechanisms remain unresolved, our findings line up with a growing body of work targeted at demystifying the complex connectivity patterns of human cerebral cortex.
Conflict of Interest: None declared.
Data for the confirmatory analysis were provided by the HCP, and the Washington University, University of Minnesota, and Oxford University Consortium (Principal Investigators David Van Essen and Kamil Ugurbil; grant 1U54MH091657) funded by 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research, and the McDonnell Center for Systems Neuroscience at Washington University. S.S.G. was partially supported by NIH grants 1R01EB020740-01A1, 1P41EB019936-01A1, 3R01MH092380-04S2, and 1U01MH108168-01. A.G. was supported by a Humboldt research fellowship from the Alexander von Humboldt Foundation.