Dissociations Between Microstructural and Functional Hierarchies within Regions of Transmodal Cortex

While the role of cortical microstructure in organising neural function is well established, it remains unclear how structural constraints can give rise to more flexible elements of cognition. While non-human primate research has demonstrated a close structure-function correspondence, the relationship between microstructure and function remains poorly understood in humans, in part because of the reliance on post mortem analyses which cannot be directly related to functional data. To overcome this barrier, we developed a novel approach to model the similarity of microstructural profiles sampled in the direction of cortical columns. Our approach was initially formulated based on an ultra-high-resolution 3D histological reconstruction of an entire human brain and then translated to myelin-sensitive MRI data in a large cohort of healthy adults. This novel method identified a system-level gradient of microstructural differentiation traversing from primary sensory to limbic regions that followed shifts in laminar differentiation and cytoarchitectural complexity. Importantly, while microstructural and functional gradients described a similar hierarchy, they became increasingly dissociated in transmodal default mode and fronto-parietal networks. Meta analytic decoding of these topographic dissociations highlighted associations with higher-level aspects of cognition such as cognitive control and social cognition. Our findings demonstrate a relative decoupling of macroscale function from cortical microstructure in transmodal regions, which likely contributes to the flexible role these regions play in human cognition.


INTRODUCTION
A core principle of neuroscience is that brain structure governs ongoing function. For decades, nonhuman primate work has confirmed the intrinsic relationship between microstructure and macrolevel function (1)(2)(3). While the role of structure in cortical processing is well characterised in the sensorymotor domain, it is less clear how this constraint gives rise to more flexible elements of cognition. In an attempt to describe how cortical areas are able to take on more abstract functional roles, Mesulam (1998) postulated a hierarchical axis of large-scale cortical organisation and connectivity, referred to as the "sensory-fugal gradient". This axis describes gradual transitions at the whole-cortex level, running from primary sensory and motor regions that are involved in externally focussed computations towards transmodal cortices that are increasingly engaged in more perceptually decoupled operations. Unlike sensory cortices, transmodal cortices have a less hierarchical organisation (3,(5)(6)(7) with a more parallel architecture that allows spatially distributed areas to respond flexibly to different types of information (8).
The current study systematically examined the interplay between microstructural constraints and functional connectivity in humans, and its contribution to high-level cognition and behaviour. Studying cortical microstructure in humans, neuroanatomists have traditionally used cell staining techniques to map spatial variations in both cyto-and myeloarchitectural features of post mortem brains (9)(10)(11)(12)(13). Extending upon this work, the 'structural model' proposes a principled relationship between microstructural similarity and large-scale connectivity of cortical areas (14), which was however primarily formulated based on non-human primate data. Notably, its holds greater predictive power in sensorimotor and unimodal association areas, compared to transmodal areas performing integrative cognitive operations. Although post mortem methods are the gold standard for describing microstructure per se, they cannot be mapped directly to function in vivo, making it hard to directly quantify how these metrics relate to neural function. With the advent of high-field magnetic resonance imaging (MRI), it has become possible to probe microstructural properties of different cortical regions in the living human brain. In particular, myelin sensitive imaging contrasts can differentiate regions with distinct myeloarchitectural profiles at an individual level (15)(16)(17)(18)(19). In parallel, resting-state functional MRI analysis can identify highly reproducible intrinsic networks formed by cortical areas (20-25). These studies have highlighted that transmodal cortex is largely composed of two spatially distributed yet functionally cohesive networks-the fronto-parietal network thought to respond to the demands of the moment in a flexible way (26-28) and the default mode network that depends on abstract information from self-generated memory and thought processes (29-31).
Core to our analysis was the formulation of a novel approach that modelled cortico-cortical networks from similarity of microstructural profiles sampled across thousands of points and in the direction of cortical columns. The model was first developed on an ultra-high resolution post mortem 3D histological reconstruction of an entire human brain (32), and we show robust evidence for a principal spatial axis of gradual cytoarchitectural variation running from primary sensory to transmodal areas, recapitulating established spatial trends in laminar differentiation and cytoarchitectural complexity (4,33). We then translated our approach to myelin-sensitive MRI in a large cohort of healthy adults, showing consistency in vivo and across individuals (34). In addition to showing correspondence between histological and MRI-derived topographies, microstructural gradients only partially converged with macroscale functional topographies derived from task-free functional connectome analysis obtained in the same subjects (25). In fact, while primary sensory regions served as a common anchor of microstructural and functional gradients, the microstructural axis depicted a progression towards limbic cortices, while its functional counterpart traversed towards default mode and fronto-parietal networks. Critically, meta-analytic decoding revealed that dissociation of function from microstructure was related to patterns of higher-order thought, such as a cognitive control or social cognition. Together our analyses on the convergence and divergence of spatial trends in microstructure and function support the hypothesis that reductions in the constraining influence of microstructure in transmodal cortex is a central underlying mechanism through which flexible cognitive functions may emerge.

Formulation of the histology-based microstructure profile covariance analysis
We modelled cortico-cortical microstructural similarity across a 100μm resolution Merker-stained 3D histological reconstruction of an entire post mortem human brain [BigBrain; https://bigbrain.loris.ca/main.php; (32) (FIGURE 1A)]. Staining intensity profiles, representing neuronal density and soma size by cortical depth, were generated along 160k surface points (henceforth, vertices) for each hemisphere (FIGURE 1B). Profile residuals, obtained after correcting intensity profile data for the y-coordinate to account for measurable shifts in intensity in anterior-toposterior direction (SI FIGURE S1), were averaged within 1012 equally sized, spatially contiguous nodes (35). Pairwise correlations of nodal intensity profiles, covaried for average intensity profile, were thresholded at 0 and positive edges were log-transformed to produce a microstructure profile covariance matrix (MPCHIST); in other words, MPCHIST captures cytoarchitectural similarity between cortical areas (see SI FIGURE S2 for distribution of values).  (SD) in residual intensity at each node are displayed on the cortex (left). Cortexwide intensity profiles were calculated by systematic intensity sampling across intracortical surfaces (rows) and nodes (columns). (C) The MPCHIST matrix depicts node-wise partial correlations in intensity profiles, controlling for the average intensity profile. Exemplary patterns of microstructural similarity from primary somatosensory (S1), anterior cingulate cortex (ACC), primary visual (V1) and the temporal pole. Seed nodes are shown in white.

Development of in vivo microstructure profile covariance analysis
The microstructure profile covariance approach was adapted to in vivo data in single individuals using T1w/T2w MRI, a ratio indexing cortical microstructure and myelin (15), shown to recapitulate sensory-fugal transitions (39). Multimodal surface-matched T1w/T2w images, pial surfaces, and white matter surfaces were obtained from the minimally-preprocessed S900 release of the Human Connectome Project (34,40). We selected a total of 219 unrelated subjects and grouped these randomly into a Discovery sample (n=110, 66 females, mean±SD age=28.8±3.8 years) and a Replication sample (n=109, 62 females, mean±SD age=28.5±3.7 years). For each individual, we systematically generated intracortical surfaces using the same equivolumetric transformation algorithm as for the histological dataset (41,42), and aggregated whole cortex intensity profiles across 64,984 linked vertices that were subsequently parcellated into 1012 contiguous nodes (35). We computed pairwise partial correlations between nodal intensity profiles (controlling for the average intensity profile), kept only positive correlations, and log-transformed the result to produce a cortexwide microstructure profile covariance matrix (MPCMRI: see SI FIGURE S5 for distribution of values). A group-average MPCMRI matrix was calculated across all participants in the Discovery sample. While microstructural similarity estimated in vivo was stronger between proximal nodes, the variance accounted for by geodesic distance was low (adjusted R 2 = 0.04, β=-2.57, p<0.001).
Diffusion map embedding revealed a principal gradient of microstructural differentiation accounting for 13.7% of variance (G1MRI, FIGURE 3B; see SI FIGURE S6 for the second gradient). In line with its histological counterpart, G1MRI was anchored on one end by primary sensory areas and on the other end by limbic regions. Cortex-wide analysis demonstrated a high correlation of G1HIST and G1MRI (r=0.63, p<0.001), driven by the close spatial correspondence of gradient extremes (FIGURE 3C, left). G1MRI depicts increasing mean myelin content, as well as a gradual transition in the relative myelin content around the midsurface (FIGURE 3B, right). Microstructural profiles of the prefrontal cortex (orange) again resembled the cortex-wide average; however, comparing node ranks across both modalities revealed a shift in prefrontal regions towards the transmodal anchor in G1MRI (FIGURE 3C, right). This effect appeared to be driven by a downward shift of lateral occipital-parietal areas towards the sensory anchor, owing to heavy myelination relative to their cytoarchitectural complexity (43). Laminar differentiation and cytoarchitectural taxonomy accounted for 44% and 37% of variance in G1MRI, respectively (SI TABLES S3-4). As in the histological analysis, the paralimbic (β=0.10, p<0.001) and idiotypic (β=0.07, p<0.001) levels were the strongest predictors within the laminar differentiation model, while motor (β=0.15, p<0.001), limbic (β=0.07, p<0.001), and primary sensory (β=0.10, p<0.001) classes were strong predictors within the cytoarchitectural model.

Correspondence of microstructural similarity with macroscale functional organisation
To examine the role of microstructural similarity in macroscale functional organisation, we generated a group-average resting state functional connectome across the Discovery subsample and derived gradients with diffusion map embedding. As shown previously (25), the principal functional gradient (G1FUNC) extends from primary sensory and motor networks, through dorsal attention and salience networks, to finally culminate in the transmodal core composed of fronto-parietal and default mode networks (FIGURE 4A). At the group level, cortex-wide analyses demonstrated moderate-to-high correlations of G1MRI with G1FUNC (r=0.52, p<0.001; FIGURE 4B) and G1HIST with G1FUNC (r=0.31, p<0.001), illustrating a common topography of microstructural profile covariance and functional connectivity. Notably, the topographies of microstructure and function were more closely related than node-to-node correspondence of microstructural similarity with functional connectivity (in vivo: r=0.10, p<0.001; histology-based: r=0.11, p<0.001), supporting the utility of connectivity-informed dimensionality reduction techniques to reveal common principles of sensory-fugal cortical organisation (44).
. All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/488700 doi: bioRxiv preprint first posted online Dec. 7, 2018; FIGURE 3. In vivo microstructure profile covariance (MPCMRI). (A) Left hemisphere pial, mid, and white matter surfaces superimposed on T1w/T2w image (left); whole-cortex intensity profiles were calculated by systematic sampling across surfaces (rows) and vertices, and then averaged with each node (columns). Mean at each node (centre right); MPCMRI matrix depicts node-wise partial correlations in intensity profiles, covaried for mean whole-cortex intensity profile (right). (B) Normalised angle matrix sorted by the principal gradient (left); variance explained by the diffusion embedding components (left centre) and the principal gradient (right centre); mean residual intensity profiles within ten discrete bins of the gradient (right). (C) Similarity of histological and in vivo gradients (G1HIST, G1MRI) shown in a density plot (left) and node-wise rank differences shown on the cortical surfaces (right). (D) Associations of G1MRI to levels of laminar differentiation (37) and cytoarchitectural class (11,38) ordered by median.
All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/488700 doi: bioRxiv preprint first posted online Dec. 7, 2018; The Replication dataset underwent identical processing procedures as the Discovery dataset. At a group level, G1MRI derived from the Replication and Discovery cohorts were highly correlated (r=0.98, p<0.001). High correlations between G1MRI and G1FUNC were also evident at the individual subject level, following alignment of individual functional and microstructural gradients to templates built from the group Discovery dataset (see METHODS, for details). For every participant, we observed a significant correlation between G1MRI and G1FUNC (0.37<r<0.61, p<0.001; FIGURE 4C-D).  (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/488700 doi: bioRxiv preprint first posted online Dec. 7, 2018; In addition to studying their commonalities, we assessed unique topographic features of the modalityspecific gradients. Cortex-wide nodal rank comparisons highlighted an upward shift in the position of the prefrontal cortex and precuneus in functional gradient space, relative to the microstructure-derived gradient, and a converse downward shift of the posterior inferior temporal and midcingulate cortices (FIGURE 5A). Notably, the greatest rank differences were evident in transmodal cortices, suggesting a specific dissociation of function from microstructure in these higher-order regions. Inspecting the distribution of nodes in intrinsic functional communities (23) along each gradient (FIGURE 5B), we noted that while the sensory-fugal gradient was overall preserved across all modalities, different sensory and transmodal networks occupied extremes in each modality. G1MRI extended from somatomotor to limbic networks, whereas G1FUNC extended from visual to transmodal default mode networks. Comparing the average node rank of each functional community, between modalities and across individuals (SI TABLE S5), we noted that the default mode and fronto-parietal networks shifted to the apex of G1FUNC, reflecting segregation of higher-order communities during rest, despite their similar myeloarchitecture. ; radar plot depicting the difference in mean node ranks of functional communities between G1MRI (blue) and G1FUNC (red), with 95% confidence intervals calculated across individuals presented with dotted lines. (C) Mean difference in the centre of gravity of meta-analysis maps in G1FUNC and G1MRI space, mapped as density plots illustrating the range of differences across subjects.
Our final analysis examined whether the functional topography diverged from microstructure specifically in cortical areas involved in abstract, perceptually decoupled functions. We conducted a meta-analysis using the Neurosynth database and estimated the centre of gravity across a set of diverse cognitive terms (25) along G1FUNC relative to G1MRI (SI FIGURE 5C). As hypothesised, the top terms exhibiting the strongest upward shift from G1MRI to G1FUNC involved multi-domain, integrative functions, such as "working memory", "verbal", "cognitive control" and "social control", whereas the strongest downward shifts involved second-order visual processing, such as "face" and "visual perception". All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

Robustness of the MPC and resultant gradients
A series of replication analyses assessed the robustness of the MPC framework. a) Histology. G1HIST was highly consistent across variations in thresholding, parcellation, intracortical surface number, and BigBrain voxel resolution (SI FIGURE S3). Furthermore, varying the α parameter for diffusion map embedding algorithm from 0 and 1 (in increments of 0.1) resulted in virtually identical gradients (all r>0.99). We also compared diffusion map derived gradients to more conventional graph theoretical characterisations of the MPCHIST matrix. Specifically, we applied Louvain modularity detection (45), incorporating a procedure that maximized consensus across fineand coarse-grained decompositions (by varying the tuning parameter γ from 0.5-1.5) across 100 repetitions. We resolved three modules (0.49<Q<0.71; SI FIGURE S7), which occupied distinct positions along the first two gradients. Diffusion map embedding thus represents a continuous and complementary approach to order microstructure profile covariance, which emphasises gradual transitions within and between discrete cortical areas that are fundamental to the hierarchical organisation of layer-specific cortical projections (46). b) In vivo. As for their histological counterpart, G1MRI solutions were robust against processing parameter variation, including matrix thresholding, parcellation scheme, and surface number (SI FIGURE S8). Again, variation of the α parameter resulted in virtually identical gradients (all r>0.99). Application of Louvain community detection identified only two modules in the microstructure profile covariance matrix which concisely halved G1MRI (0.35<Q<0.41, SI FIGURE S9). We replicated G1MRI in an independent dataset of healthy late adolescents/young adults (47) that underwent magnetisation transfer imaging (spatial correlation: r=0.78, p<0.001, SI FIGURE S10) and in a cohort of 17 healthy adults scanned at our imaging centre in whom quantitative T1 relaxation data was available (spatial correlation: r=0.81, p<0.001, SI FIGURE S11), demonstrating robustness of the approach to acquisition site, acquisition type, and surface construction.

DISCUSSION
Cortical areas identified by classical neuroanatomical studies represent discrete regions embedded within gradual transitions in cyto-and myeloarchitecture (9)(10)(11)48). Cortical gradients are nearly ubiquitous across microstructural and functional domains of the mammalian neocortex, where mounting evidence supports a common and overarching "sensory-fugal" organisation (46,49). Leveraging a ultra-high-resolution 3D reconstruction of an entire post mortem human brain (32), we sampled microstructure profiles and utilised unsupervised techniques to identify smooth transitions in cytoarchitectural composition. We translated our approach from histology to myelin-sensitive in vivo MRI and recovered similar microstructural gradients, which were consistent across individuals. Collectively, our findings support a common sensory-fugal axis of inter-areal cortical differentiation across histology and in vivo microstructure. Importantly, we observed a clear dissociation between microstructural and functional organisation that was most prominent in transmodal cortex within the default mode and fronto-parietal networks. While transmodal regions were distributed across the microstructural hierarchy, they were localised to the top of the functional hierarchy. This dissociation in transmodal regions was related to aspects of human cognition such as cognitive control and social cognition providing important support for views that these regions are relatively untethered by cortical structure and so allowing them to take on more flexible cognitive roles.
To bridge different scales of human brain organisation (50), we developed a cortico-cortical network model based on microstructural similarity between areas. Microstructure profile generation, based on state-of-the-art equivolumetric surface construction techniques (41), provided a simplified recapitulation of cellular changes across the putative laminar structure of the cortex. Importantly, our covariance framework generated cytoarchitecturally grounded networks, which were sensitive to laminar thickness as well as cell size and density. Complementing classic histological investigations that describe and delineate specific cortical areas (9-11) and more recent work identifying specific histological features of cortical areas from the BigBrain dataset (51), the current study revealed gradual cytoarchitectural transitions across the folded neocortex. Consistent with long standing views of cortical organisation from studies of non-human primates, our analysis elucidated a macroscopic All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
By showing a close correlation between individualised microstructural gradients and corresponding functional connectome gradients (25), our findings expand the structural model of connectivity, initially formulated in non-human primates (14), to living humans. Microstructural and functional gradients were consistently anchored by primary sensory regions on the lower end, recapitulating their consistent location at the bottom of functional and structural hierarchies in primates (3,4,25,33,49,55). In primary regions, where neural processes are strongly constrained by sensorimotor contingencies, information processing is driven by extrinsic inputs and interactions with the outside world (56,57), as well as intrinsic signalling molecules involved in axonal guidance, cell adhesion, regional circuit formation, and thus cortical patterning (58). The impact of these forces decreases with synaptic distance, producing a sensory-fugal axis of neurostructural and functional differentiation along the cortical mantle. In fact, given the default mode network's location as the most distant functional system relative to the location of primary sensory sulci (25), it may be relatively un-tethered by the influence of intrinsic signalling molecules and extrinsic activity (59), and thus expresses less unique microstructural profiles, as observed here and in classical neuroanatomical studies (9)(10)(11). Furthermore, in contrast to a close functional coupling of microstructurally similar areas at the low end of the cortical gradient, the functional connectivity between higher-order regions may not be as strongly constrained by microstructural similarity. For example, prefrontal and parietal association regions of the higher-order networks densely interact functionally (21,23,60,61) despite distinct cytoarchitectural compositions (62). A plausible mechanism for reduced influence of microstructural homophily in higher-order networks may be the increased relevance of long-range structural connections (63). Such wider ranging structural connections may accommodate the more distributed spatial layout and diverse functional roles of transmodal cortices. Conversely, primary sensory and motor regions exhibit more locally clustered short range connectivity profiles related to microstructural similarity (64,65), likely in accordance with their more specialized and stable functional roles.
Our approach was robust with respect to variations in algorithmic and analytical choices. While the singular nature of the BigBrain dataset prohibited replication of the histological pipeline, we demonstrated consistency of neuroimaging-derived microstructural gradients across three cohorts with unique myelin sensitive contrasts (15,17,19,34,47). Importantly, the in vivo approach can serve as a lower resolution, yet biologically meaningful extension of the histological work. In fact, myelin-based gradients exhibited a comparable association with laminar differentiation and cytoarchitectural complexity as the histology-based gradient. With the emergence of algorithms that detect cortical laminae based on histological data (51) and increasing availability of ultra-high field MRI scanners, we expect that the histological as well as in vivo MPC approach may be further refined in future work, providing an even more direct bridge between the micro-and macro-scales of human brain organisation. In addition, as it offers a hierarchy-dependent reference frame to examine the interplay of microstructural and system-level network mechanisms in single individuals, the proposed framework may be advantageous in neurodevelopmental research and complement existing covariance network mapping approaches that have shown promise in studying typical and atypical development (66)(67)(68)(69)(70)(71). As such, it offers a powerful and freely available (http://github.com/micaopen/MPC) tool to investigate coordinated changes in cortical microstructure paralleling the emergence and maturation of large-scale functional networks during development (72,73), and may conversely provide insight into the structural underpinnings of atypical network configurations in complex neurodevelopmental disorders (74).
We close by considering the significance of the observed progressive dissociation between structural and functional topographies in cortical areas involved in multi-domain, integrative processing. Components of the fronto-parietal network have been shown to guide behaviour in an adaptive All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/488700 doi: bioRxiv preprint first posted online Dec. 7, 2018; manner, changing in line with external demands (26), and possibly substantiated by its ability to dynamically reconfigure its functional connectivity (75). Although the functional role of the default mode network as the putative apex of the cortical functional hierarchy in primates remains subject to debate (76), it is evidently involved a broad class of memory-driven operations, involving selfreferential and simulative thought processes and some degree of abstraction (29-31). The interesting possibility that reduced microstructural constraints enable functional diversity and flexibility was also supported by ad hoc meta-analytical synthesis, which suggested that structure-function decoupling preferentially related to working memory, social cognition, and cognitive control. Such a hypothesis provides a potential mechanistic account for why some of the more creative acts of the human mind emerge through the interaction of the two most dominant yet structurally diverse functional systems (77). All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

Histological data acquisition and preprocessing
An ultra-high resolution Merker stained 3D volumetric histological reconstruction of a post mortem human brain from a 65-year-old male was obtained from the open-access BigBrain repository on February 2, 2018 [https://bigbrain.loris.ca/main.php; (32)]. The post mortem brain was paraffinembedded, coronally sliced into 7,400 20μm sections, silver-stained for cell bodies (78) and digitised. Manual inspection for artefacts (i.e., rips, tears, shears, and stain crystallisation) was followed by automatic repair procedures, involving non-linear alignment to a post mortem MRI, intensity normalisation, and block averaging (79). 3D reconstruction was implemented with a successive coarse-to-fine hierarchical procedure (80). We downloaded the 3D volume at four resolutions, with 100, 200, 300, and 400μm isovoxel size. We primarily analysed 100μm data and used 200, 300 and 400μm data to assess consistency of findings across spatial scales. Computations were performed on inverted images, where staining intensity reflects cellular density and soma size. Geometric meshes approximating the outer and inner cortical interface (i.e., the GM/CSF boundary and the GM/WM boundary) with 163,842 matched vertices per hemisphere were also available (81).

Histology-based microstructure profile covariance (MPC) analysis
a) Surface sampling. We systematically constructed 10-100 equivolumetric surfaces, in steps of 1, between the outer and inner cortical surfaces (42). The equivolumetric model compensates for cortical folding by varying the Euclidean distance, ρ, between pairs of intracortical surfaces throughout the cortex to preserve the fractional volume between surfaces (41). ρ was calculated as follows for each surface where α represents fraction of the total volume of the segment accounted for by the surface, while Aout and Ain represent the surface area of the outer and inner cortical surfaces, respectively. Next, vertexwise microstructure profiles were estimated by sampling intensities along linked vertices from the outer to the inner surface across the whole cortex. In line with previous work (51), layer 1 was approximated as the top ten percent of surfaces and removed from the analysis due to little interregional variability. Note however, that findings were nevertheless virtually identical when keeping the top ten percent of surfaces. To reduce the impact of partial volume effects, the deepest surface was also removed. Surface-based linear models, implemented via SurfStat for Matlab (http://www.math.mcgill.ca/keith/surfstat/) (82), were used to account for an anterior-posterior increase in intensity values across the BigBrain due to coronal slicing and reconstruction (32), whereby standardised residuals from a simple linear model of surface-wide intensity values predicted by the midsurface y-coordinate were used in further analyses.
b) MPC matrix construction. Cortical vertices were parcellated into 1,012 spatially contiguous cortical 'nodes' of approximately 1.5cm 2 surface area, excluding outlier vertices with median intensities more than three scaled median absolute deviations away from the node median intensity. The parcellation scheme preserves the boundaries of the Desikan Killany atlas (35) and was transformed from conte69 surface to the BigBrain midsurface via nearest neighbour interpolation. Nodal intensity profiles underwent pairwise Pearson product-moment correlations, controlling for the average whole-cortex intensity profile. MPCHIST for a given pair of nodes i and j was thus where rij is the Pearson product moment correlation coefficient of the BigBrain intensity profiles at nodes i and j, ric the correlation coefficient of the intensity profile at node i with the average intensity profile across the entire cortex, and rjc the Pearson correlation of the intensity profile at node j with the average intensity profile across the whole brain. The MPC matrix was absolutely thresholded at All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/488700 doi: bioRxiv preprint first posted online Dec. 7, 2018; zero then remaining MPC values were log-transformed to produce a symmetric 1012☓1012 MPCHIST matrix. The in-house developed code for MPC construction is available online (https://github.com/MICA-MNI/micaopen/MPC). c) Parameter estimation. The optimal surface number was determined based on the stability of the MPC matrix. This procedure involved (repeatedly and randomly) dividing the vertex intensity profiles within each node into two groups and constructing two MPC matrices, then calculating the Euclidean distance between them. The procedure was repeated 1000 times. Although the MPC matrix instability was robust to variations in surface number, the 18-surface solution exhibited a noticeable local minimum MPC instability in the studied range (10-100 surfaces) and was used in subsequent analyses (SI FIGURE S2). Notably, the MPC gradient was similar using two finer grained solutions (i.e., 54 and 91 surfaces), where local minima were observed as well. More details on the origins of the stability statistic in clustering algorithms may be found elsewhere (83). d) Relation to spatial proximity. To determine whether MPCHIST was not purely driven by spatial proximity, we correlated MPCHIST strength with geodesic distance for all node pairs. The latter was calculated using the Fast Marching Toolbox between all pairs of vertices, then averaged by node (https://github.com/gpeyre/matlab-toolboxes/tree/master/).

Histology-based MPC gradient mapping
In line with previous studies (25,84), the MPCHIST matrix was proportionally thresholded at 90% per row, then converted into a normalised angle matrix. We then applied diffusion map embedding (36), a nonlinear manifold learning technique, to identify principal gradient components explaining MPCHIST variance in descending order (each of 11012). In brief, the algorithm estimates a low-dimensional embedding from a high-dimensional affinity matrix. In this space, cortical nodes that are strongly inter-connected by either many suprathreshold edges or few very strong edges are closer together, whereas nodes with little or no inter-covariance are farther apart. The name of this approach, which belongs to the family of graph Laplacians, derives from the equivalence of the Euclidean distance between points in the embedded space and the diffusion distance between probability distributions centred at those points. Compared to other non-linear manifold learning techniques, the diffusion maps algorithm is relatively robust to noise and computationally inexpensive (85,86). Notably, the algorithm is controlled by a single parameter α, which controls the influence of the density of sampling points on the manifold (α = 0, maximal influence; α = 1, no influence). In this and previous studies (25,84), we followed previous recommendations and set α=0.5, a choice that retains the global relations between data points in the embedded space and has been suggested to be relatively robust to noise in the covariance matrix. Gradients were mapped back onto BigBrain midsurface and visualised with SurfStat for Matlab (http://www.math.mcgill.ca/keith/surfstat/) (82), and we assessed the amount of MPCHIST variance explained. To show how the principal gradient in MPCHIST (G1HIST) relates to systematic variations in microstructure, we calculated and plotted the mean microstructure profiles within ten equally-sized discrete bins of G1HIST.

Relation of G1HIST to laminar differentiation and cytoarchitectural taxonomy
We evaluated correspondence of G1HIST to atlas information on laminar differentiation and cytoarchitectural class. To this end, each cortical node was assigned to one of four levels of laminar differentiation (i.e., idiotypic, unimodal, heteromodal or paralimbic) derived from a seminal model of Mesulam, which was built on the integration of neuroanatomical, electrophysiological, and behavioural studies in human and non-human primates (37), and one of the seven Von-Economo/Koskinas cytoarchitectural classes (i.e., primary sensory, secondary sensory, motor, association 1, association 2, limbic or insular) (11,38). In the case of laminar differentiation maps, assignment was done manually; in the case of cytoarchitectural classes, we mapped previously published Von Economo/Koskinas classes (71) to the BigBrain midsurface with nearest neighbour interpolation and assigned nodes to the cytoarchitectural class most often represented by the underlying vertices. Finally, we estimated the contribution of level of laminar differentiation (D, a All rights reserved. No reuse allowed without permission. (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
Structural and resting state functional MRI data underwent HCP's minimal preprocessing (40,87,88). For structural MRI, images underwent gradient nonlinearity correction. When repeated scans were available, these were co-registered and averaged. Following brain extraction and readout distortion correction, T1w and T2w images were co-registered using rigid body transformations. Subsequently, non-uniformity correction using T1w and T2w contrast was applied (89). Preprocessed images were nonlinearly registered to MNI152 space and cortical surfaces were extracted using FreeSurfer 5.3.0-HCP (90-92), with minor modifications to incorporate both T1w and T2w (15). Cortical surfaces in individual subjects were aligned using MSMAll (93,94) to the hemisphere-matched conte69 template (95). T1w images were divided by aligned T2w images to produce a single volumetric T1w/T2w image per subject (Glasser and Van Essen, 2011). Notably, this contrast nullifies inhomogeneities related to receiver coils and increases sensitivity to intracortical myelin.
For resting state functional MRI, the timeseries were corrected for gradient nonlinearity and head motion. The R-L/L-R blipped scan pairs were used to correct for geometric distortions. Distortion corrected images were warped to T1w space using a combination of rigid body and boundary-based registrations (96). These transformations were concatenated with the transformation from native T1w to MNI152, to warp functional images to MNI152. Further processing removed the bias field (as calculated for the structural image), extracted the brain, and normalised whole brain intensity. A highpass filter (>2000s FWHM) corrected the time series for scanner drifts, and additional noise was removed using ICA-FIX (97). Tissue-specific signal regression was not performed (98,99). Following the minimal preprocessing pipeline, we transformed resting state functional MRI to native space and sampled time-series at each vertex of the MSMAll-registered (93,94) mid-thickness cortical surfaces.

In vivo microstructure profile covariance (MPC) analysis
We estimated MPC in the in vivo dataset (MPCMRI) in the same manner as MPCHIST with the only adjustment that intensity profiles were not corrected for y-coordinates; instead, the contrast reversal for T1w and T2w data was used to correct for inhomogeneity as part of the minimal processing pipeline of the Human Connectome Project. We generated equivolumetric surfaces between the outer and inner cortical surfaces (see Equation (1)), and systematically sampled T1w/T2w values along 64,984 linked vertices from the outer to the inner surface across the whole cortex. In turn, MPCMRI can be denoted as an extension of Equation (2), in which MPCMRI(i,j) for a given pair of nodes i and j is defined by: where s is a subject and n is the number of subjects. We systematically evaluated matrix stability across different intracortical surface numbers and selected the most stable solution (SI FIGURE S8).

In vivo MPCMRI gradient: relation to the sensory-fugal gradient and the histological MPC gradient
As for the histological data, diffusion map embedding derived a principal gradient (G1MRI) from the group average MPCMRI matrix. Correspondence of the in vivo G1MRI to the histological G1HIST gradient was estimated via Spearman rank correlation between spatially matched nodes. To localise differences, we calculated the difference in rank of each node. As before, we assessed the contribution of level of laminar differentiation and cytoarchitectural class to the G1MRI via multiple regression.

Correspondence of microstructure and functional connectivity
Individual functional connectomes were generated by averaging preprocessed timeseries within nodes, correlating nodal timeseries and converting them to z scores. For each individual, the four available resting state scans were averaged at the matrix level. Then, a group average functional connectome was calculated across the Discovery cohort. Correlation coefficients were calculated between the group average functional connectome and the MPCHIST and MPCMRI matrices. The group average functional connectome was proportionally thresholded at 90% per row, transformed into a cosine similarity matrix, transformed into a normalised angle matrix, then diffusion map embedding was applied, producing G1FUNC. We calculated the correspondence of the G1FUNC with G1HIST and G1MRI with Spearman rank correlations. The differences in the gradients were localised by comparing node ranks across the whole cortex and within functional communities. Seven functional communities were mapped onto the conte69 surfaces from a previously published parcellation (23) with nearest neighbour interpolation from fsaverage5, then nodes were assigned to the functional community most often represented by the underlying vertices. To aid interpretation of the modality-specific gradients, G1MRI and G1FUNC were discretised into fifty equally sized bins and we calculated the proportion of each bin accounted for by each functional community, then performed seven paired t-tests contrasting average node rank of a functional community in G1MRI and G1FUNC across individuals.
We used the NeuroSynth meta-analytic database (www.neurosynth.org) (100) to assess which cognitive functions were associated with shifts between G1MRI and G1FUNC. To this end, we downloaded and surface projected meta-analytic z-statistic maps of 24 terms covering a wide range of cognitive functions, which correspond to the topic names defined by Margulies et al., (2016). We discretised G1MRI and G1FUNC into five-percentile bins and, for each term, calculated the mean zstatistic within each bin. From this, we deduced the centre of gravity of each term within gradient space, and contrasted the term-specific centre of gravity in G1MRI and G1FUNC at an individual level.

Individualised MPCMRI gradients and relation to the individual-specific functional hierarchy
Inter-individual consistency of G1MRI was assessed in the Replication dataset. The MPCMRI pipeline was deployed at an individual level, thus resolving individualised MPCMRI matrices and gradients. Additionally, the diffusion map embedding was employed on functional connectomes to derive individual functional gradients (25). To ensure the spatial correspondence of individual gradients, the individual gradients from the Replication dataset underwent Procrustes linear alignment to the All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/488700 doi: bioRxiv preprint first posted online Dec. 7, 2018; Discovery dataset group average embedding. Individual cross-modal coupling was calculated as the Spearman rank correlation between G1MRI and G1FUNC.

Robustness of G1HIST
We assessed the robustness of G1HIST by altering pipeline parameters, repeating MPCHIST generation and diffusion map embedding, then calculating the Pearson correlation of the modified G1HIST with the original gradient. In particular, we evaluated variable matrix thresholds (i.e., 70-95%, in steps of 1%), alternative surface number in which MPCHIST matrix instability reached a local minima (i.e., 54and 91-surface solutions), the voxel resolution of the BigBrain volume (100-400 um), and spatial scale (i.e., vertex-vs parcel wise construction); for the latter, we correlated nodal gradient values with the median vertex within each parcel.

Robustness of G1MRI
We repeated the robustness procedures reported for the histological gradients in the in vivo dataset, including variation of thresholding level, parcellation usage and surface number. Here, the pipeline was repeated with 23 surfaces, pertaining to local minima in the in vivo MPCMRI matrix instability.
Cortical surfaces were extracted from the T1-weighted scans using FreeSurfer 6.0 (90-92), and fourteen equivolumetric intracortical surfaces were generated (42). qT1 was registered to Freesurfer native space using a boundary-based registration (96), and Freesurfer native space was registered to standard conte69 space using Caret5 landmark-based registration (95). We used the former to sample qT1 intensity values along the intracortical surfaces, and the latter to resample the evaluated surfaces to a common conte69 space with 64,984 matched vertices. As in the main approach, we averaged vertex-wise intensity profiles within 1012 nodes (35), computed pairwise partial correlations between nodal intensity profiles (controlling for the average intensity profile), kept only positive correlations, and log-transformed the result to produce a MPCMRI-QT1 matrix. Finally, we generated a group-average MPCMRI-QT1 matrix and applied diffusion map embedding. The similarity of G1MRI-QT1 to the original G1MRI was measured with a node-wise Spearman rank correlation.  (72)). We similarly applied the MPC framework to the MT profiles, then averaged the MPCMRI-MT matrix across the group and applied diffusion map embedding. We measured the Spearman rank correlation between G1MRI-MT and the original G1MRI, which was recalculated with 308 matched cortical areas.
All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. FIGURE S1-11 TABLE S1-5 All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/488700 doi: bioRxiv preprint first posted online Dec. 7, 2018; FIGURE S1. Horizontal view of the volumetric reconstruction of the "BigBrain" with pial, mid and white matter surfaces projected on the left hemisphere. Notably, we corrected for the linear relationship between intensity values and midsurface y-coordinate (r=-0.68, p<0.001), which existed due to coronal slicing and reconstruction of the BigBrain.
All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/488700 doi: bioRxiv preprint first posted online Dec. 7, 2018; BigBrain datasets were characterised as lower resolution replications, and G1HIST was found to be highly correlated across these resolutions (all r>0.81, all p<0.001).
All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

FIGURE S4. (A)
The second principal component, accounting for 12.7% of variance in MPCHIST components, is projected on the BigBrain midsurface. (B) Scatterplot depicting the first two embedding gradients, with corresponding probability density functions. The second gradient divides the lower-order areas of the first gradient, insomuch that somatomotor and primary visual areas (red) are separated from ventral prefrontal areas and secondary visual areas (blue).
All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

FIGURE S6. (A)
The second principal component, accounting for 11.7% of variance in MPCMRI components, projected on the conte69 midsurface. (B) Scatterplot depicting the first two embedding gradients, with corresponding probability density functions. The second gradient divides the higher-order areas of the first gradient, insomuch that the cingulate, orbitofrontal cortex and the inferior temporal gyrus (red) are separated from the prefrontal cortex, precuneus, temporo-parietal junction and superior temporal gyrus (blue).
All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/488700 doi: bioRxiv preprint first posted online Dec. 7, 2018; (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/488700 doi: bioRxiv preprint first posted online Dec. 7, 2018; 1 Frontal and temporal association areas, displayed in yellow (FIGURE 2) 1 Parietal and superior temporal association, displayed in purple (FIGURE 2) All rights reserved. No reuse allowed without permission.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
(which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.