A multi-modal parcellation of human cerebral cortex (original) (raw)

Nature. Author manuscript; available in PMC 2017 Feb 11.

Published in final edited form as:

PMCID: PMC4990127

EMSID: EMS68870

Matthew F Glasser,1 Timothy S Coalson,#1 Emma C Robinson,#2,3 Carl D Hacker,#4 John Harwell,1 Essa Yacoub,5 Kamil Ugurbil,5 Jesper Andersson,2 Christian F Beckmann,6,7 Mark Jenkinson,2 Stephen M Smith,2 and David C Van Essen1

Matthew F Glasser

1Department of Neuroscience, Washington University Medical School, Saint Louis, Missouri 63110, USA

Timothy S Coalson

1Department of Neuroscience, Washington University Medical School, Saint Louis, Missouri 63110, USA

Emma C Robinson

2FMRIB centre, Nuffield Department of Clinical Neurosciences, John Radcliffe Hospital, University of Oxford, Oxford OX3 9DU, UK

3Department of Computing, Imperial College, London SW7 2AZ, UK

Carl D Hacker

4Department of Biomedical Engineering, Washington University, Saint Louis, Missouri 63110, USA

John Harwell

1Department of Neuroscience, Washington University Medical School, Saint Louis, Missouri 63110, USA

Essa Yacoub

5Center for Magnetic Resonance Research (CMRR), University of Minnesota, Minneapolis, Minnesota 55455, USA

Kamil Ugurbil

5Center for Magnetic Resonance Research (CMRR), University of Minnesota, Minneapolis, Minnesota 55455, USA

Jesper Andersson

2FMRIB centre, Nuffield Department of Clinical Neurosciences, John Radcliffe Hospital, University of Oxford, Oxford OX3 9DU, UK

Christian F Beckmann

6Donders Institute for Brain, Cognition and Behavior, Radboud University, Nijmegen 6525 EN, The Netherlands

7Department of Cognitive Neuroscience, Radboud University Medical Centre Nijmegen, Postbus 9101, Nijmegen 6500 HB, The Netherlands

Mark Jenkinson

2FMRIB centre, Nuffield Department of Clinical Neurosciences, John Radcliffe Hospital, University of Oxford, Oxford OX3 9DU, UK

Stephen M Smith

2FMRIB centre, Nuffield Department of Clinical Neurosciences, John Radcliffe Hospital, University of Oxford, Oxford OX3 9DU, UK

David C Van Essen

1Department of Neuroscience, Washington University Medical School, Saint Louis, Missouri 63110, USA

1Department of Neuroscience, Washington University Medical School, Saint Louis, Missouri 63110, USA

2FMRIB centre, Nuffield Department of Clinical Neurosciences, John Radcliffe Hospital, University of Oxford, Oxford OX3 9DU, UK

3Department of Computing, Imperial College, London SW7 2AZ, UK

4Department of Biomedical Engineering, Washington University, Saint Louis, Missouri 63110, USA

5Center for Magnetic Resonance Research (CMRR), University of Minnesota, Minneapolis, Minnesota 55455, USA

6Donders Institute for Brain, Cognition and Behavior, Radboud University, Nijmegen 6525 EN, The Netherlands

7Department of Cognitive Neuroscience, Radboud University Medical Centre Nijmegen, Postbus 9101, Nijmegen 6500 HB, The Netherlands

#Contributed equally.

Supplementary Materials

Extended Supplementary Methods.

GUID: 156B8075-5AB4-427B-B6CD-3004ABFFD60F

Neuroanatomical Supplementary Results.

GUID: 32DFD9A7-088D-465B-9BC7-7E8E2255F18F

Supplementary Information Guide.

GUID: 73852F79-81C6-47ED-B151-9C879137D563

Supplementary Results and Discussion.

GUID: 9DCC9C6D-DC5D-481F-968D-E6E4F77A21F3

Abstract

Understanding the amazingly complex human cerebral cortex requires a map (or parcellation) of its major subdivisions, known as cortical areas. Making an accurate areal map has been a century-old objective in neuroscience. Using multi-modal magnetic resonance images from the Human Connectome Project (HCP) and an objective semi-automated neuroanatomical approach, we delineated 180 areas per hemisphere bounded by sharp changes in cortical architecture, function, connectivity, and/or topography in a precisely aligned group average of 210 healthy young adults. We characterized 97 new areas and 83 areas previously reported using post-mortem microscopy or other specialized study-specific approaches. To enable automated delineation and identification of these areas in new HCP subjects and in future studies, we trained a machine-learning classifier to recognize the multi-modal ‘fingerprint’ of each cortical area. This classifier detected the presence of 96.6% of the cortical areas in new subjects, replicated the group parcellation, and could correctly locate areas in individuals with atypical parcellations. The freely available parcellation and classifier will enable substantially improved neuroanatomical precision for studies of the structural and functional organization of human cerebral cortex and its variation across individuals and in development, aging, and disease.

Neuroscientists have long sought to subdivide the human brain into a mosaic of anatomically and functionally distinct, spatially contiguous areas (cortical areas and subcortical nuclei), as a prerequisite for understanding how the brain works. Areas differ from their neighbours in microstructural architecture, functional specialization, connectivity with other areas, and/or orderly intra-area topographic organization (for example, the map of visual space in visual cortical areas)1–3. Accurate parcellation provides a map of where we are in the brain, enabling efficient comparison of results across studies and communication among investigators; as a foundation for illuminating the functional and structural organization of the brain; and as a means to reduce data complexity while improving statistical sensitivity and power for many neuroimaging studies.

The human cerebral cortex has been estimated to contain anywhere from ~50 (ref. 1) to ~200 (refs 3, 4) areas per hemisphere. However, attaining a consensus whole-cortex parcellation has been difficult because of practical and technical challenges that we address here.

Most previous parcellations were based on only one neurobiological property (such as architecture, function, connectivity or topography), and many cover only part of the cortex. Combining multiple properties provides complementary as well as confirmatory information, as different properties distinguish different sets of areal boundaries, and more confidence can be placed in boundaries that are consistent across multiple independent properties. We analysed all four properties across all of neocortex in both hemispheres, using new or refined methods applied to the uniquely rich repository of exceptionally high-quality magnetic resonance imaging (MRI) data provided by the Human Connectome Project (HCP), which benefited from major advances in image acquisition and preprocessing5–8. Architectural measures of relative cortical myelin content and cortical thickness were derived from T1-weighted (T1w) and T2-weighted (T2w) structural images5,9,10. Cortical function was measured using task functional MRI (tfMRI) contrasts from seven tasks11. Resting-state functional MRI (rfMRI) revealed functional connectivity of entire cortical areas plus topographic organization within some areas.

Previous parcellations typically used either fully automated algorithmic approaches, or else manual or partly automated neuroanatomical approaches in which neuroanatomists delineated areal borders, documented areal properties, and identified areas after consulting prior literature. Here we combined both approaches. For the initial parcellation, we adapted a successful observer-independent semi-automated neuroanatomical approach for generating post-mortem architectonic parcellations12,13 to non-invasive neuroimaging. We used an algorithm to delineate potential areal borders (transitions in two or more of the cortical properties described above), which two neuroanatomists (authors M.F.G. and D.C.V.E.) then interpreted, documenting areal properties and identifying areas relative to the extant neuroanatomical literature. We then used a fully automated algorithmic approach, training a machine-learning classifier to delineate and identify cortical areas in individual subjects based on multi-modal areal fingerprints, allowing the parcellation to be replicated in new subjects and studies.

Prior parcellations have either used small numbers of individuals or group averages that are ‘blurry’ from inaccurate alignment of brain areas across subjects. We aligned cortical data using ‘areal features’, including maps of relative myelin content and resting state networks that are more closely tied to cortical areas than are the folding patterns typically used for alignment14. The markedly improved intersubject cortical alignment using cortical folding, myelin and resting state fMRI enabled us to generate the ‘typical subject’s’ parcellation from a highly detailed 210-subject group average data set.

Group map reproducibility of fine details

We analysed two independent groups of HCP subjects—210P (‘parcellation’) and 210V (‘validation’)—aligned using areal-feature-based registration (called ‘MSMAll’, see Methods section on image preprocessing). Figure 1 illustrates the consistency of fine spatial patterns for maps reflecting relative myelin content (left panels) and a task-fMRI language-related activation (right panels). The maps are strikingly similar across the 210P and 210V groups, including variations in relative myelin content within the primary somatosensory cortex related to somatotopic organization (white and black arrows in left panels, see legend) and small features in the task fMRI maps (white ellipses in right panels). Supplementary Results and Discussion 1.1 with associated Supplementary Figs 1–5 include more examples of such cross-group consistency for architecture (myelin and cortical thickness), function (tfMRI contrast maps), and two resting state connectivity measures.

An external file that holds a picture, illustration, etc. Object name is emss-68870-f001.jpg

Consistency of fine spatial details in independent group averages.

Relative myelin content maps (left hemisphere) and task fMRI contrast beta maps from the LANGUAGE story contrast (right hemisphere) on inflated (columns 1 and 3) and flattened surfaces (columns 2 and 4). Rows 1 and 2 are the group averages of the 210P and 210V data sets, respectively. White and black arrows indicate consistent variations in myelin content within primary somatosensory cortex that are correlated with somatotopy (see Supplementary Neuroanatomical Results 6 and Supplementary Neuroanatomical Results Fig. 8). The white oval indicates a small, sharp, and reproducible feature in the right hemisphere of the LANGUAGE story contrast. Relative myelin content will hereafter be referred to as myelin (see legend of Supplementary Fig. 1 in Supplementary Results and Discussion 1.2). Data at (http://balsa.wustl.edu/WDpX).

Because of areal-feature-based alignment, maps of average cortical folding in this study are much blurrier than are maps of areal properties because folding patterns and area locations are imperfectly correlated4,12,14 (for example, compare Supplementary Fig. 7e and 7j (group average and individual subject folding) in Supplementary Results and Discussion 1.3). Group average folding patterns remain sharp mainly in early sensory areas, where areal locations and folds are tightly correlated (for example, central and calcarine sulci, see Supplementary Fig. 1, rows 3 and 4, in Supplementary Results and Discussion 1.1). The regional difference in sharpness between maps of areal properties and folding highlights the importance of alignment based on areal features, rather than folding patterns, as a prerequisite for accurately parcellating group average data. The high spatial resolution of the HCP’s MRI images and lack of aggressive spatial smoothing10 prior to group averaging also contribute to making our maps substantially sharper than those from traditional neuroimaging studies.

Quantitatively, the 210P and 210V group average datasets were highly correlated across the cortical surface (r = 0.998 for myelin, r = 0.994 for cortical thickness after correction for folding-related effects; r = 0.996 and r = 0.979 for two folding-related measures (FreeSurfer’s ‘sulc’ and ‘curv’, respectively); r = 0.995, r = 0.984 and r = 0.944 for the maximum, median and minimum, respectively, of the task fMRI contrasts, and a median reproducibility of r = 0.989 for two measures of resting state connectivity). These excellent map reproducibilities provide confidence that the parcellation will reflect the areal pattern of typical subjects in the healthy young adult population. See Methods section on modalities for parcellation, and Supplementary Methods 1.3–3.4 for the methods used to generate these maps.

A 180-area group average parcellation

To identify transitions representing candidate areal boundaries, we designed and implemented a semi-automated, quantitative approach adapted for multi-modal neuroimaging data represented on two-dimensional cortical surface models (see Methods section on the gradient-based parcellation approach and Supplementary Methods 4.1–5.3). The approach is similar in spirit to a highly successful semi-automated observer-independent approach13,15. However, instead of objectively identifying potential areal borders in postmortem histological sections, we identified them algorithmically on the cortical surface by computing the first derivative of each areal feature map (its spatial gradient magnitude)16. Candidate borders were then interpreted by the neuroanatomists to exclude artefacts. Each area’s properties were documented (in the Supplementary Neuroanatomical Results), and putative areas were related to the extant neuroanatomical literature.

These semi-automated approaches contrast with classical observer-dependent parcellation approaches1,3 that relied on visual inspection to locate often subtle transitions in cortical architecture and with some modern observer-dependent retinotopic parcellation methods17,18. They also differ from fully automated, unsupervised methods19–21 in which the outcomes depend heavily on algorithmic input parameters (for example, thresholds or number of requested clusters) and are not validated by a neuroanatomist.

Area 55b illustrates our multi-modal gradient-based parcellation approach using gradients of three areal feature maps (see Fig. 2). Area 55b is a small, elongated, and notably distinct area (outlined in black or white) bounded by the frontal eye field (FEF) and premotor eye field (PEF), primary motor cortex (4), ventral premotor cortex (6v), and prefrontal areas 8Av and 8C. In the myelin map (Fig. 2a), area 55b is lightly myelinated and lies between moderately myelinated areas FEF (above) and PEF (below), just anterior to heavily myelinated primary motor cortex (area 4). Thus, area 55b is surrounded on three sides by myelin gradients (Fig. 2e). Area 55b is strongly activated in the ‘Story versus Baseline’ task contrast from the HCP’s ‘LANGUAGE’ task (Fig. 2b) and is entirely surrounded by a strong gradient for this task contrast (Fig. 2f). It also has distinctive functional connectivity, as revealed by a seed location (lightly myelinated area PSL) selectively connected with 55b (Fig. 2c) and a different seed location (heavily myelinated area LIPv) strongly connected with FEF and PEF (Fig. 2d) but not with 55b. The result is strong mean gradients in dense functional connectivity surrounding 55b (Fig. 2g). Ref. 22 illustrated area 55b on a schematic surface map (Fig. 2h) as a lightly myelinated area bounded on three sides by more heavily myelinated areas. Because of the similarity to the dorsal portion of 55b in ref. 22, we use the same name.

An external file that holds a picture, illustration, etc. Object name is emss-68870-f002.jpg

Parcellation of exemplar area 55b using multi-modal information.

The border of 55b is indicated by a white or black outline. a, Myelin map. b, Group average beta map from the LANGUAGE Story versus Baseline task contrast. c, d, Functional connectivity correlation maps from a seed in area PSL (white sphere, arrow) (c) and a seed in area LIPv (white sphere, arrow) (d). e, Gradient magnitude of the myelin map shown in a. f, Gradient magnitude of the LANGUAGE Story versus Baseline task contrast shown in b. g, Mean gradient magnitude of the functional connectivity dense connectome (see section on modalities for parcellation in the Methods). h, A dorsal schematic view of the prefrontal cortex as parcellated in ref. 22, in which shading indicates the amount of myelin found using histological stains of cortical grey matter. Data at (http://balsa.wustl.edu/Qv4P).

To generate the complete parcellation of 180 areas and area complexes in each hemisphere, we adopted a systematic, objective, and quantitative approach (see the gradient-based parcellation approach section in the Methods and in Supplementary Methods 5.1–5.3). Our major criteria, met in nearly all cases, included: (i) spatially overlapping gradient ‘ridges’ between each pair of areas for at least two independent areal feature maps; (ii) similar gradient ridges present in roughly corresponding locations in both hemispheres; (iii) gradients that were not correlated with artefacts; and (iv) robust and statistically significant cross-border differences in the feature maps. Another consideration (but not a requirement) was whether published evidence exists for a boundary in an approximately corresponding location. Studies with publicly available parcellations registered onto atlas surfaces4 were directly compared with our data; however, most regions required indirect comparisons with published figures (for example, Fig. 1h). Initial areal boundaries meeting these criteria were delineated by two neuroanatomists (authors M.F.G. and D.C.V.E.).

In a second computational stage, the path of each manually drawn border was optimized algorithmically using gradients of the most informative feature maps selected by the neuroanatomists (those with visually obvious gradients and differences across the border). These feature maps were confirmed to have robust and statistically significant differences across the final border. The semi-automated gradient-based parcellation approach is further described in Supplementary Methods 5.1–5.3), and the entire semi-automated process is illustrated for area V1 in Supplementary Neuroanatomical Results 1; other sections of this document describe and illustrate the information used to delineate and the literature used to name all 180 cortical areas.

Figure 3 shows the multi-modal cortical parcellation in the left and right hemispheres on inflated and flattened surfaces, with areal boundaries delineated by black contours. A total of 180 areas and areal complexes per hemisphere is near the higher end of earlier estimates noted above3,4. We consider 180 likely to be a lower bound, as some parcels are probably complexes of multiple areas (for example, based on finer-grained published parcellations, and other regions suffer from reduced sensitivity due to fMRI signal loss). Many areas (83/180) were assigned names based on published parcellations from dozens of separate studies that used a variety of invasive or specialized methods (see Table 2 in Supplementary Neuroanatomical Results), reflecting how far the field has been from a consensus neuroanatomical parcellation. Some of the newly described 97 areas have de novo names (for example, DVT for the dorsal visual transitional area), while others represent finer-grained parcellations of previously reported areas (for example, area 31 into areas 31a, 31pd, and 31pv). A few represent complexes in which a published finer grained parcellation was not visible in our data (for example, areas 29 and 30 combined into area the retrosplenial complex (RSC)), but these may be again subdivided once higher resolution data is available. The 180 areas differ widely in their shapes, sizes, and the positions of their borders relative to cortical folds.

An external file that holds a picture, illustration, etc. Object name is emss-68870-f003.jpg

The HCP’s multi-modal parcellation, version 1.0 (HCP_MMP1.0).

The 180 areas delineated and identified in both left and right hemispheres are displayed on inflated and flattened cortical surfaces. Black outlines indicate areal borders. Colours indicate the extent to which the areas are associated in the resting state with auditory (red), somatosensory (green), visual (blue), task positive (towards white), or task negative (towards black) groups of areas (see Supplementary Methods 5.4). The legend on the bottom right illustrates the 3D colour space used in the figure. Data at http://balsa.wustl.edu/WN56).

The parcellation in Fig. 3 is coloured to reflect each area’s degree of association in the resting state (determined using multiple regression, see Supplementary Methods 5.4) with five functionally specialized groups of areas: early auditory (red), early somatosensory/motor (green), and early visual areas (blue). These represent the three dominant input streams to the brain. Also used were two core groups of cognitive areas that are strongly anti-correlated in our data, the task positive network (towards white) and task negative (also called the default mode) network (towards black). Hence, the strongly bluish, greenish, and reddish regions are predominantly but not exclusively associated with visual, somatosensory-motor, and auditory processing, respectively. Qualitatively, the predominantly unimodal regions appear to collectively occupy less than half of the neocortical sheet. Areas probably more strongly biomodal include blue-green areas such as LIPv and MT (visual and somatosensory-motor) and purple areas such as POS2 and RSC (visual and auditory). The remaining regions form a complex mosaic, with some intermixing of lighter (task-positive) and darker (task-negative) areas along with many light or dark pastel hues suggestive of ‘cognitive’ areas that may be preferentially associated with one or another sensory modality. The bilateral symmetry of functional organization is striking, in that nearly all areas have qualitatively very similar hues in the left and right hemispheres. However, interesting colour asymmetries occur in a few areas, especially language-related areas 55b, PSL, SFL, and 44 and their right hemisphere homologues, which also have asymmetric task-fMRI functional profiles (see Supplementary Neuroanatomical Results 8, 15, 21 and 22).

Internal heterogeneity is evident in some cortical areas, particularly those with topographically organized representations. In the somatosensory-motor strip (largely architecturally defined somatosensory and motor areas 3a, 3b, 1, 2, and 4), we identified five clearly defined topographic subareas in resting state and task fMRI data (see Supplementary Neuroanatomical Results 6 and the associated Supplementary Fig. 8). In this parcellation we treat topographic subdivisions as ‘subareas’ rather than calling them full ‘areas’. For visual cortex, its visuotopic organization revealed a set of hemifield representations in each hemisphere, something not achieved in previous unsupervised resting state functional connectivity-based parcellations19–21,23. Also, ultrahigh-field MRI reveals sub-areal cortical organization along both laminar24,25, and columnar26,27 axes, so our parcellation represents one of many important levels of granularity in brain organization.

Cross-validation of the parcellation

The initial statistical analysis used in the semi-automated parcellation was circular, to the extent that the 210P dataset was used for both creating and testing the parcellation. Hence, we carried out an additional statistical cross-validation using the 210V dataset and a comprehensive set of feature maps (see the statistical cross validation of the multi-modal parcellation section of the Methods and Supplementary Results and Discussion 1.2). This analysis also reveals which areal properties were most useful in defining areal boundaries (a condensed representation of the detailed information provided in the Supplementary Neuroanatomical Results). Supplementary Fig. 6 in Supplementary Results and Discussion 1.2 shows four independent categories of features: cortical thickness, myelin maps, task fMRI, and resting state fMRI, and how many of these categories showed robust and statistically significant differences across each areal border. Fully 96% of areal borders had robust effect sizes (Cohen’s d > 1) in two or more feature categories and all were statistically significant after correcting for multiple comparisons in two or more feature categories in cross-border, across-subject _t_-tests. Resting state fMRI was the most useful category, followed by task fMRI, myelin maps, and lastly cortical thickness, which was consistent with the neuroanatomists’ observations and documentation in the Supplementary Neuroanatomical Results.

Exemplar parcellation-based analyses

Spatial smoothing is often used to increase the signal-to-noise ratio (SNR) in neuroimaging analyses, to try to compensate for inaccurate registration of brain areas, and/or to satisfy statistical assumptions. However, smoothing blurs data across boundaries between areas (on the surface) and tissue compartments (in the volume). An areal parcellation enables area-wise analyses (averaging data within each area), thereby improving SNR and statistical power without the deleterious effects of spatial smoothing (to the extent that properties within an area are uniform). Parcellation dramatically reduces data dimensionality, illustrated here using the HCP’s myelin, thickness, task, and resting state data (Fig. 4).

An external file that holds a picture, illustration, etc. Object name is emss-68870-f004.jpg

Example parcellated analyses using the HCP’s multi-modal cortical parcellation.

a, Dense myelin maps on lateral (top) and medial (bottom) views of inflated left hemisphere. b, c, Example dense (b) and parcellated (c) task fMRI analysis (LANGUAGE story versus baseline) expressed as Z statistic values. d, The entire HCP task fMRI battery’s Z statistics for 86 contrasts (47 unique, see section on modalities for parcellation in the Methods) analysed in parcellated form and displayed as a matrix (rows are parcels, columns are contrasts, white outline indicates the map in c). e, A major improvement in Z statistics from fitting task designs on parcellated time series instead of fitting them on dense time series and then parcellating afterwards (blue points are 360 parcels × 86 task contrasts; note the upward tilting deviation from the red line). f, Parcellated myelin maps. g, A parcellated folding-corrected cortical thickness map (in mm). h, i, Parcellated functional connectivity maps on the brain (seeded from area PGi, black dot). These parcellated connectomes are computed using either full or partial correlation (see Supplementary Methods 7.1). In both cases, the task negative (default mode) network is apparent. j, A parcellated connectome matrix view with the full correlation connectome below and the partial correlation connectome above the diagonal (white line shows the displayed partial correlation brain map). Data at (http://balsa.wustl.edu/RG0x).

The ‘dense’ (vertex-wise) myelin map shown in Fig. 4a has ~30,000 surface grey matter vertices per hemisphere, whereas a ‘parcellated’ myelin map (Fig. 4f) shows the same overall pattern with 180 cortical areas (vertices within an area have the same value, see also Fig. 4g for parcellated cortical thickness). Example dense and parcellated task fMRI analysis contrast maps (Figs. 4b, c, LANGUAGE Story versus Baseline) can be represented as a single column (white) in a 180-area by 86-task-contrast matrix (Fig. 4d). Parcellated analyses hold great promise for task fMRI studies, as they improve the signal-to-noise ratio by averaging fMRI time series within parcels prior to fitting the task design, increasing Z statistics (Fig. 4e). Parcellation is effectively a neurobiologically constrained smoothing approach that also increases statistical power by efficiently consolidating otherwise non-independent statistical tests. This approach will benefit studies aimed at understanding the functional and structural organization of the brain in health or disease at an area-wise level (studies that currently summarize results using three-dimensional coordinates in a standardized stereotaxic space). Parcellated analyses also aid in the clarity and efficiency of communicating results (for example: “area 55b in the left hemisphere showed a statistically significant +1% BOLD activation in my language task”).

Parcellated analyses are comparably useful when characterizing structural or functional connectivity, as previously recognized23,28. Preprocessing of HCP data results in fMRI data represented as ‘grayordinates’ (cortical grey matter surface vertices and subcortical grey matter voxels5). A dense connectome, containing connectivity between all pairs of 91,282 grayordinates is ~3.3 × 105-fold larger than an area-wise parcellated connectome for ~500 areas (connectivity between all pairs of areas), yet the parcellated connectome captures the neurobiologically relevant variance at the areal level. Parcellated connectomes are illustrated using a seed location in area PGi (black dot) for full correlation (Fig. 4h) and partial correlation (Fig. 4i) functional connectivity brain maps together with their associated parcellated connectome matrices (Fig 4j, full correlation below and partial correlation above the diagonal). In both cases, the task negative (default mode) network is evident, though the partial correlation connectome is much sparser than the full correlation connectome.

Individuals with atypical areal patterns

The precisely aligned group average multi-modal cortical parcellation represents the overall spatial arrangement of cortical areas in the ‘typical’ individual from a healthy young adult population. However, we found atypical topological arrangements of some areas in some individuals that are discernible across multiple modalities, including resting-state networks, task-fMRI activations, and myelin maps. Distinguishing genuinely atypical areal topologies from inadequately aligned typical patterns depended on the MSMAll areal-feature-based registration to align cortical areas precisely. We summarize key findings here and extensively characterize this important phenomenon in the Supplementary Results and Discussion 1.3.

Previously described area 55b and neighbouring areas FEF and PEF showed particularly notable individual differences in topological arrangements. For the 210P subjects, 89% showed the typical configuration in the left hemisphere (area 55b bordered by area PEF inferiorly and area FEF superiorly, as in Fig. 2), which was well aligned with group average area 55b after MSMAll registration. However, in one subgroup (4%, n = 9), a patch having the multi-modal characteristics of area 55b is shifted superiorly relative to the upper limb subregion of sensori-motor cortex (Supplementary Figure 7 in Supplementary Results and Discussion 1.3). In another subgroup (6%, n = 12), area 55b is split into two pieces by a merger of areas FEF and PEF, rather than the typical splitting of FEF and PEF by 55b (Supplementary Fig. 8 in Supplementary Results and Discussion 1.3). Such topological deviations in individual subjects’ areal maps raise intriguing questions for future exploration. They also cannot be corrected by a topology-preserving registration aimed at aligning individual subjects’ areas with the group average ‘atlas’ parcellation. Thus, we introduce an alternative fully automated cortical parcellation approach that can identify and delineate both typical and atypical areas in individual subjects that were not a part of the original 210P group.

Automated individual-subject parcellation

The semi-automated neuroanatomical approach described above is impractical for de novo individual subject parcellation of all ~1,100 HCP subjects having complete MRI datasets so as to identify the atypical areal topologies mentioned above. Instead, we developed an automated method for generating individual subject parcellations based on a supervised machine learning classifier previously used to identify resting state functional networks in individual subjects29. In our case, the areal classifier learns the multi-modal ‘areal fingerprint’ of each cortical area that distinguishes it from surrounding cortex. Based on multi-modal feature maps that represent the areal properties of architecture, function, connectivity, and topography, the areal classifier returns a prediction (0% to 100%) that each area exists at a given cortical surface vertex. The highest prediction value across areas at each vertex is used to generate the individual subject parcellation (see the cortical areal classifier section in Methods and in Supplementary Methods 6.1–6.8). Once trained using the 210P subjects (and a separate ‘29T’ group of test subjects, see the subjects and acquisitions section of the Methods), the areal classifier should be able to use only the multi-modal areal fingerprints that it has learned to reproduce the parcellation in an independent group of validation subjects (210V).

A critical early test of the areal classifier was whether it could accurately and reliably map areas that are not aligned with the population-based atlas parcellation after MSMAll areal-feature-based alignment (see Supplementary Results and Discussion 1.4). Examples of successful classification of areas 55b, FEF, and PEF are shown in Supplementary Fig. 9 of the Supplementary Results and Discussion for typical subjects, shifted 55b subjects, and split 55b subjects. In each illustrated case, the classifier correctly identified 55b and its neighbours (as assessed by the neuroanatomists’ inspection of the multi-modal areal features shown in the figure). Supplementary Fig. 10 in the Supplementary Results and Discussion 1.4 shows that these atypical 55b topologies and classifications are stable across widely spaced repeat scanning sessions in a ‘test–retest’ group of 27 subjects (see Methods section on subjects and acquisition).

Areal detection and parcellation consistency

Another critical test of both the parcellation and the areal classifier is the classifier’s performance in detecting the 180 cortical areas in individual subjects, particularly in independent validation subjects that were not used to generate the parcellation or train the classifier. The top two rows of Fig. 5 show the performance of the classifier in detecting each area (see the cortical areal classifier section in Methods). Importantly, the classifier aims to detect whole areas based on their multi-modal fingerprints, rather than detecting differences in areal features across paired areal boundaries as was done in the cross-validation analysis (Supplementary Fig. 6 in Supplementary Results and Discussion 1.2). The overall areal detection rate was 98.0% of all areas across all subjects for the 210P parcellation and training dataset (row 1) and 96.6% for the independent 210V validation dataset (row 2), indicating excellent overall performance of the areal classifier.

An external file that holds a picture, illustration, etc. Object name is emss-68870-f005.jpg

Areal detection rates, probabilistic areas, and parcellation reproducibility.

Rows 1 (210P) and 2 (210V) show the individual subject areal detection rates (see Methods section on cortical areal classifier) as parcellated maps. Most areas are yellow (100%), and the minimum detection rate across both rows was 73%. Rows 3 and 4 illustrate probabilistic maps of areas V1, M1, RSC, MT, LIPv, TE1a, 46, and 10r for the 210P (row 3) and 210V (row 4) groups. Row 5 shows the original parcellation derived from the semi-automated neuroanatomical approach. Row 6 shows the group MPM maps from 210P (blue), 210V (red), and their overlap (purple). Data at (http://balsa.wustl.edu/WL8m).

The areal classifier was used to generate probabilistic maps of each cortical area (illustrating residual variability in spatial location after MSMAll areal feature-based registration), and to assess the reproducibility of the parcellation in the independent 210V dataset. Rows 3 and 4 of Fig. 5 show strikingly similar probabilistic maps of 8 non-overlapping areas with differing degrees of spatial variability (V1, M1, RSC, MT, LIPv, TE1a, 46, and 10r) from the 210P and 210V groups. All probability maps were combined to produce a group maximum probability map (MPM), where the area with the highest probability at each vertex was found. Row 5 shows the original semi-automated parcellation borders, and row 6 compares the group MPM maps from 210P (blue) and 210V (red), with purple representing overlapping vertices. The borders in Row 6 are almost entirely purple, indicating very high reproducibility of the group MPM maps (r = 0.965, dice = 0.960, see the cortical areal classifier section of Methods). This reproducibility is similar to that of the original group average feature maps discussed above. The correlation of the original semi-automated parcellation (row 5) with the 210P group MPM (row 6) was r = 0.913, dice = 0.902, indicating that the classifier made modest adjustments to better fit the data. We predict there will be very high reproducibility of the parcellation across the rest of the ~1,100 subject HCP dataset. Example individual subject parcellations and their reproducibility based on repeated scan sessions are shown in Supplementary Fig. 11 of Supplementary Results and Discussion 1.4. The individual parcellations are reasonably reproducible (median r = 0.77, dice = 0.72) but, unsurprisingly, not as reproducible as the group parcellations, which benefit from averaging across many subjects. Other analyses yield interesting information about the sizes of cortical areas in the group average and variability in areal size across individuals (Supplementary Results and Discussion 1.5).

Generalizing the classifier for future studies

In contrast to the semi-automated approach (above) where neuroanatomists chose the information to delineate and identify the 180 cortical areas in group-average data (see Supplementary Neuroanatomical Results), the areal classifier automatically determines (without human intervention) what information is most useful for delineating and identifying these cortical areas in individual subjects. As illustrated in Supplementary Figure 12 of Supplementary Results and Discussion 1.6, the areal classifier uses the task fMRI data least, perhaps because task fMRI feature maps are noisier in individual subjects than other feature maps and its information content is largely redundant with the resting state data30. This finding is important for generalizability of the areal classifier to other studies because replicating the customized, hour-long HCP task fMRI battery is unlikely to be feasible for most neuroimaging studies. Ideally, the areal classifier would be able to perform nearly as well relying only on architecture, connectivity, and topography. Accordingly, we trained the classifier again on the 210P dataset, but omitted the task fMRI-based feature maps. When trained this way, the classifier indeed performed nearly as well as when all features were used, detecting 97.6% of areas in 210P (versus 98.0% using all features) and 96.4% of areas in 210V (versus 96.6% using all features). Hence, we anticipate that the areal classifier will generalize to other studies that acquire the following core set of MRI images: high-resolution T1w and T2w; spin echo-based b0 field map; and extensive fMRI data acquired using ‘multiband’ pulse sequences to improve spatial and temporal resolution7 (see Supplementary Results and Discussion 2.3 and 2.9). These are the same image acquisition requirements as the HCP’s minimal preprocessing pipelines5 and the MSMAll areal feature-based registration pipeline14 (Supplementary Methods 2.4). Future studies adhering to these image acquisition guidelines will be able to use the unified framework of the HCP’s analysis pipelines to automatically generate individualized parcellated analyses from unprocessed MRI images, a major advance over traditional neuroimaging methods that have often relied on comparisons with Brodmann’s hand drawn parcellation published in 1909 (ref. 1).

Discussion

We have produced a population-based 180-area per hemisphere human cortical parcellation using exceptionally high quality multimodal data from hundreds of Human Connectome Project subjects aligned using an improved areal feature-based cross-subject alignment method (MSMAll). Inspired by an observer-independent post-mortem architectural parcellation approach13, we developed a semi-automated neuroanatomical approach adapted to non-invasively acquired multi-modal MRI data. Though algorithms determined the final areal borders, the multi-modal data were carefully interpreted by neuroanatomists, the properties of each cortical area were documented, and each area was named in relation to the extant neuroanatomical literature (see Supplementary Neuroanatomical Results). A cross-validation showed that the areas forming the parcellation were robustly and statistically significantly different from their neighbours across multiple modalities. We identify this parcellation as HCP-MMP1.0 (Human Connectome Project Multi-Modal Parcellation version 1.0), with version 1.0 anticipating future refinements as better data become available (see Supplementary Results and Discussion 2.1).

Unexpectedly, we discovered that despite improved intersubject alignment, some areas have atypical topological arrangements in some subjects, which we demonstrated for areas 55b, FEF, and PEF. We developed a fully automated method for parcellating individual subjects based on a machine learning classifier that can cope with this kind of individual variability. The areal classifier detected 96.6% of individual subject cortical areas in new subjects, including atypical areas, and replicated the group parcellation in an independent sample. Though we made extensive use of the HCP’s specialized task fMRI battery when generating the parcellation, we showed that task fMRI data is not essential for future studies aiming to use the areal classifier to automatically define the cortical areas in their subjects. Instead, it suffices to acquire the same core set of MRI images needed for the rest of the HCP’s software pipelines.

By generating a robust neuroanatomical map of human neocortical areas—a century-old aim of neuroscience—and providing methods for mapping these areas in any individual undergoing study with non-invasive neuroimaging, the present work stands apart from previous human cortical parcellations. The overall approach described here shows that we can produce sharp, reproducible brain images across multiple non-invasive neuroimaging modalities. We can generate a highly reproducible and generalizable cortical parcellation through state-of-the-art methods of data acquisition, preprocessing, and analysis designed to compensate for individual variability and thereby minimize blurring of images. These improvements, together with the new parcellation, make it desirable to use spatial localization methods that move beyond the traditional use of stereotaxic coordinates combined with Brodmann areal assignments to characterize centers of cortical activation in fMRI studies. From a neuroanatomical perspective, there has often been substantial uncertainty whether any two neuroimaging studies have found results in the same cortical areas or not. The situation is analogous to astronomy in which ground-based telescopes produced relatively blurry images of the sky before the advent of adaptive optics and space telescopes.

Many topics are discussed further in the Supplementary Results and Discussion 2.1–2.10 (for example, avenues for improving the parcellation and other issues left for future work, further discussion of the neuroscientific implications of our results, and additional datasets that could profitably be linked to our parcellation). As the topographic organization of higher cognitive areas becomes better understood, some parcels currently considered to be full areas may later be considered to be subareas of larger topographically organized cortical areas (analogous to somatotopic subregions of topographically organized sensory and motor areas illustrated in Supplementary Neuroanatomical Results 6). Though our use of multiple modalities probably mitigates this issue relative to traditional uni-modal parcellations, the extent to which the human multi-modal cortical parcellation may be revised along such lines remains a question for future work using the state of the art methods mentioned above (see Supplementary Results and Discussion 2.8).

The MSMAll registration and the areal classifier are or will soon be freely available on GitHub, the visualization tool Connectome Workbench (http://humanconnectome.org) and the parcellation, data, and scenes for reproducing each of the figures are in the BALSA database31. These tools can be used to address fundamental questions regarding the identification of cortical areas, when reporting results or thinking about and discussing brain organization in relation to studies of human cognition, lifespan, and disease. Several additional interesting avenues of investigation are now open. The ability to discriminate individual differences in the location, size, and topology of cortical areas from differences in their activity or connectivity should facilitate the dissection of how each property is related to behaviour and genetic underpinnings, for example, in learning disabilities or those with distinctive cognitive traits. The ability to non-invasively and automatically delineate cortical areas in living subjects may have clinical implications, for example by providing neurosurgeons with detailed, individualized maps of the brains on which they operate. There are also important implications for our understanding of human cortical evolution. The dramatic expansion in neocortex along the human lineage occurred mainly in higher cognitive regions of lateral prefrontal, parietal, and temporal cortices9,13,32,33. Comparisons with nonhuman primates, including marmosets and macaques (both widely used in invasive studies), and great apes, may yield new insights regarding the emergence of new cortical areas and the divergences in areal functions, which collectively led to the cognitive capabilities that make us uniquely human as a species and as individuals.

Methods

Subjects and acquisition

A total of 449 young adult twins and non-twin siblings (ages 22–35) from the Human Connectome Project (HCP) were scanned according to the HCP’s acquisition protocol5–7. The MRI acquisition included collecting T1w and T2w structural images, task-based and resting state-based fMRI images, diffusion-weighted images, and b0 field maps. Images were acquired at high spatial and temporal resolution on a customized Siemens 3 tesla (3T) scanner and with customized slice accelerated sequences for fMRI (see Supplementary Methods 1.1–1.2). All subjects from the HCP 500-subject data release (July, 2014) having complete fMRI sessions were included. They were divided into two independent groups of 210 subjects that shared no family members between them, together with a remaining group of 29 test (29T) subjects that shared family members with 210P but not 210V. The first group of subjects (210P, 130 females, 80 males) was used for creating the parcellation and training the areal classifier, which also made use of the 29T group to avoid overfitting. The second group of subjects (210V, 116 females, 94 males) was used only for statistical cross-validation of the parcellation, areal classifier detection rates in independent subjects, and group parcellation reproducibility measures. A test–retest group of 27 subjects scanned twice through the entire MRI protocol and independently processed through the HCP pipelines was used for individual subject reproducibility measures. Subject recruitment procedures and informed consent forms, including consent to share de-identified data, were approved by the Washington University institutional review board. Datasets were de-identified and are publicly shared on the ConnectomeDB database (https://db.humanconnectome.org).

Image preprocessing

Spatial image preprocessing (distortion correction and image alignment) was carried out using the HCP’s spatial minimal preprocessing pipelines5. This included steps to maximize alignment across image modalities, to minimize distortions relative to the subject’s anatomical space, and to minimize spatial smoothing (blurring) of the data. The data were projected into the 2 mm standard CIFTI grayordinates space, which includes cortical grey matter surface vertices and subcortical grey matter voxels5. This offers substantial improvements in spatial localization over traditional volume-based analyses, enabling more accurate cross-subject and cross-study registrations and avoiding smoothing that mixes signals across differing tissue types or between nearby cortical folds. Additionally, we did minimal smoothing within the CIFTI grayordinates space to avoid mixing across areal borders prior to parcellation.

For cross-subject registration of the cerebral cortex, we used a two-stage process based on the multimodal surface matching (MSM) algorithm14 (see Supplementary Methods 2.1–2.5). An initial ‘gentle’ stage, constrained only by cortical folding patterns (FreeSurfer’s ‘sulc’ measure), was used to obtain approximate geographic alignment without overfitting the registration to folding patterns, which are not strongly correlated with cortical areas in many regions. Previously, we found that more aggressive folding-based registration (either MSM-based or FreeSurfer-based) slightly decreased cross-subject task-fMRI statistics, suggesting that aligning cortical folds too tightly actually reduces alignment of cortical areas14. A second, more aggressive stage used cortical areal features to bring areas into better alignment across subjects while avoiding neurobiologically implausible distortions or overfitting to noise in the data. The areal features used were myelin maps, resting state network maps computed with weighed regression (an improvement over dual regression34 described in the Supplementary Methods 2.3) and resting state visuotopic maps (see Supplementary Methods 4.4). Areal distortion was measured by taking the log base-2 of the ratio of the registered spherical surface tile areas to the original spherical surface tile areas. The mean (across space) of the absolute value of the areal distortion averaged across subjects from both registration stages was 30% less than the standard FreeSurfer folding-based registration and the maximum (across space) of this measure was 54% less. Despite less overall distortion, the areal-feature-based registration delivers substantially more accurate registration of cortical areas than does FreeSurfer folding-based registration as judged by cross-subject task fMRI statistics, an areal feature that was not used to drive the registration14. Because MSM registration preserves topology and is relatively gentle (it does not tear or distort the cortical surface in neurobiologically implausible ways), it is unable to align some cortical areas in some subjects where the areal arrangement differs from the group average (see Supplementary Results and Discussion 1.3–1.4 for more details on atypical areas). Group average registration drift away from the gentle folding-based geographic alignment was removed from the surface registration35 (see Supplementary Methods 2.5) to enable comparisons of this dataset with datasets registered using different areal features (for example, post-mortem cytoarchitecture). Group average registration drift is any consistent effect of the registration during template generation on the mean size, shape, or position of areas on the sphere (as opposed to the desired reductions in cross-subject variation). An obvious example is the 37% increase in average brain volume produced by registration to MNI space4. Uncorrected drifts during surface template generation can cause apparent changes in cortical areal size, shape, and position when comparing across studies.

Resting state fMRI data were denoised for spatially specific temporal artefacts (for example, subject movement, cardiac pulsation, and scanner artefacts) using the ICA+FIX approach, which includes detrending the data and aggressively regressing out 24 movement parameters36,37. We avoided regressing out the ‘global signal’ (mean grey-matter time course) from our data because preliminary analyses showed that this step shifted putative connectivity-based areal boundaries so that they lined up less well with other modalities, likely because of the strong areal specificity of the residual global signal after ICA+FIX clean up. Task fMRI data were temporally filtered using a high pass filter. More details on resting state and task fMRI temporal preprocessing are described in the Supplementary Methods 1.6–1.8. Substantial spatial smoothing was avoided for both datasets, and all images were intensity normalized to account for the receive coil sensitivity field. Artefact maps of large vein effects, fMRI gradient echo signal loss, and surface curvature were computed as described in Supplementary Methods 1.9.

Modalities for parcellation

The multi-modal cortical parcellation used information related to the four areal properties of architecture, function, connectivity, and topography2. Architecture was measured using T1w/T2w myelin content maps plus cortical thickness maps with surface curvature regressed out5,9,10 (Supplementary Methods 1.5). Function was measured using task-fMRI responses to 7 tasks in 86 task contrasts (47 unique; 39 were sign-reversed contrasts). Effect size maps (beta maps) after correction for the receive field were used instead of Z statistic maps because we are interested in regional differences in the magnitude of the BOLD (blood oxygen level dependent) signal change induced by the tasks, rather than differences in the significance of the BOLD signal change. Functional connectivity was measured using pairwise Pearson correlation of the denoised resting state time series of each pair of grayordinates. Topographic organization was explored using resting state time series in visual cortex, with spatial regressors representing polar angle and eccentricity patterns in area V1 combined with a modified ‘dual-regression-like’ approach that weights each surface vertex according to the cortical surface area that it represents (see Supplementary Methods 4.4). The semi-automated multi-modal parcellation was generated using group average data for all of these modalities from the 210P group of subjects (see Supplementary Methods 3.1–3.3 for details on how the group averages were created for each modality). The reproducibility of these group average maps was assessed by correlating the spatial maps for the 210P and 210V groups (see Supplementary Results and Discussion 1.1).

The gradient-based parcellation approach

Classically, cortical areas have been defined based on sharp changes in one or more of the areal properties of architecture, function, connectivity, and topography. Traditionally, this relied heavily on visual inspection, until more objective and quantitative approaches became available12,13. One highly successful approach to post-mortem architectural parcellation involves computing a dissimilarity metric, (the Mahalanobis distance) between neighbouring feature profiles generated from segmented histological images and testing for statistically significant and large spikes in dissimilarity that indicate putative areal boundaries. For in vivo data, a similarly powerful approach involves taking the first derivative (the spatial gradient) of a measure of interest along cortical surface and using the gradient magnitude to objectively identify locations where the measure is changing rapidly. One can then draw putative areal boundaries along the resulting gradient ridges10,16,19. Here we combined elements of both approaches in a multi-modal context to generate semi-automatically drawn areal borders that were then evaluated statistically. Gradients were computed for architectural, functional, connectivity, and topographic modalities (see Supplementary Methods 4.1–4.4).

To incorporate expert knowledge and priors from the neuroanatomical literature into the parcellation process, the neuroanatomists (authors M.F.G and D.C.V.E.) evaluated the multi-modal neuroimaging data and its gradients to define initial areal borders based on the following criteria. (1) Presence of a co-localized gradient ridge in at least two independent modalities was taken as strong evidence of an areal border, and the vast majority of areal borders satisfied this criterion. (2) Presence of corresponding gradients in the left and right hemispheres provided further evidence for a genuine areal border. For the vast majority of borders, the same modalities yielded robust gradients in both hemispheres. We did not find strong evidence for an area present in one hemisphere that was absent in the other (though a few areas show hemispheric asymmetries in their functional ‘signature’ and/or in their spatial relationships with neighbouring areas). (3) We ignored gradients clearly attributable to imaging artefacts (see Supplementary Neuroanatomical Results for details). (4) Cortex on opposite sides of the border needed to differ robustly and significantly in the areal features used to delineate the border. (5) Confidence was increased if prior literature described a corresponding areal border. (6) Early runs of a supervised machine learning algorithm (see the cortical areal classifier section of the Methods below) needed to be able to learn to distinguish each cortical area from its neighbours in a large majority of individual subjects based on individual subject multi-modal features (the early runs were only done using 210P and 29T, keeping 210V independent for later analyses). After the neuroanatomists delineated the initial areal borders and chose the important areal features that defined them, an automated algorithm then optimized the border placement so that it followed the most probable path based on the chosen areal features (see Supplementary Methods 5.1–5.3). The Supplementary Neuroanatomical Results documents the information that was used to distinguish each of the 180 areas from its neighbours.

The neuroanatomists named areas based on previous parcellations whenever a reasonable match to the literature could be made. In some cases, areal identification was based on the similarity of the area’s properties relative to previously reported areas (for example, area 4, primary motor cortex, is known to be heavily myelinated and thick; area V2 has a mirror-image visuotopic map relative to neighbouring area V1). In most cases, however, the information used to describe previous cortical areas (for example, cytoarchitecture) was not available in the HCP data, and areal identification mainly reflected spatial correspondences relative to cortical folding patterns (if reliable for that region of cortex) or spatial relationships between neighbouring cortical areas. The strongest evidence for areal identification came from studies that provided surface-based probabilistic or maximum probability maps, ideally also registered using areal features and dedrifting of templates35. In these cases, we directly compared these data with our data and show the degree of overlap in the Supplementary Neuroanatomical Results. When such data were unavailable, we used published information to the degree feasible (see Supplementary Methods 5.3 for limitations of non-surface-based/not publicly available data) to make areal identifications or to describe new areas that had not previously been identified. The information used to name each cortical area is described in the Supplementary Neuroanatomical Results.

Statistical cross validation of the multi-modal parcellation

Once the parcellation has been created, parcellated representations of data from each modality can be generated using either the group parcellation or the individual subject parcellations. For the statistical cross-validation, we created parcellated myelin, cortical thickness, task fMRI, and resting state functional connectivity datasets using the semi-automated multimodal group parcellation (see Supplementary Methods 7.1). For myelin and cortical thickness, we simply averaged the values of the dense individual subject maps within each area. For task fMRI, we averaged the time series within each area prior to computing task statistics (to benefit from the CNR improvements of parcellation demonstrated in Fig. 4e). For the same reason, we averaged resting state time series within each parcel prior to computing functional connectivity to form a parcellated functional connectome.

For each pair of areas that shared a border in the parcellation, we computed a paired samples two-tailed _t_-test across subjects on these parcellated data for each feature (ignoring tests that involved the diagonal in the resting state parcellated functional connectome). We thresholded these tests at the Bonferroni-corrected significance level of P < 9 × 10−8 (number of area pairs across both hemispheres (1,050) × number of features (266) × number of tails (2) × 0.05) and an effect size threshold of Cohen’s _d_ > 1. We grouped the features into 4 independent categories (cortical thickness, myelin, task fMRI, and resting state fMRI) to determine for each area pair whether it showed robust and statistically significant differences across multiple modalities. For more details, see Supplementary Methods 7.2.

The cortical areal classifier

We used a supervised machine learning classifier to automatically delineate and identify each cortical area from its neighbours across a large majority of individual subjects based on multi-modal information. Besides validating the robustness of the parcellation, this provides useful information about each individual subject’s parcellation, along with an approach to generalizing the parcellation to other datasets. To automatically parcellate individual subjects, we adapted the multi-layer perceptron used by ref. 29 to delineate and identify seven resting state networks more accurately than simpler linear methods including dual regression. We used the multi-layer perceptron to classify all 180 areas in our parcellation using multi-modal feature maps and relied on two neuroanatomically sensible assumptions to simplify the problem. (1) After areal feature-based registration (MSMAll), we assumed that each cortical area was approximately in the same general location across subjects (for example, we don’t expect to find V1 outside the occipital lobe). This also means that we consider widely separated regions having similar multi-modal areal fingerprints to be distinct cortical areas even if they have similar architecture, coactivation in functional tasks, and belong to the same resting state network. These assumptions allowed us to reduce the overall classification problem to a set of 180 classification problems per hemisphere, each involving discrimination of one area from the areas around it. (2) Also, instead of classifying each area from all of its neighbours specifically (one class for the area plus one class for each neighbouring area), we set up the problem as a binary classification (the most robust kind of classification problem), classifying each area from all of the surrounding cortex as a single alternate class. This surrounding cortex represents a ‘searchlight’ for the area, and this searchlight was the group parcel location plus a 30 mm radius surrounding the group parcel in all directions across the surface (meaning that for a 10 mm circular area, the searchlight would be a circle of 70 mm in diameter, still a quite large region of cortex). The 30 mm radius (geodesic distance computed on the group average mid-thickness surface corrected for vertex area loss due to averaging) was chosen because it easily encompassed the individual variation in area 55b in the 210P group (55b approaches a worst case because it is a relatively small and highly variable cortical area). The training labels were the group area from the semi-automated parcellation (class 1), and the remaining cortex in the searchlight (class 2).

The features used by the classifier covered the same set of modalities used for the original parcellation, including architectural measures of myelin and cortical thickness with curvature regressed out; task fMRI maps (redundant information was reduced and Supplementary Neuroanatomical Results increased with a d = 20 ICA-based approach run on the task contrast beta maps, see Supplementary Methods 6.4); the 77 surface-related resting state fMRI network maps computed on individual subjects using weighted regression from an overall d = 137 group ICA; five visuotopic topographic maps transformed into a format interpretable by the classifier; and maps of artefacts that the classifier used to interpret differences in areal features due to artefactual effects (see Supplementary Methods 6.3–6.5 for further description of each modality’s classifier features). These 112 multi-modal feature maps were generated for each vertex in each of the 449 subjects and the 27 repeated subjects, with each hemisphere processed separately. Other than the 30 mm radius searchlight region of interest (ROI), the classifier has no spatial concept of where the area should be (it operates independently on each vertex and only knows what the area’s fingerprint looks like in the feature space). Consequently, special consideration was given to the spatial visuotopic patterns, which were transformed into maps whose values reflected the alternating mirror symmetric organization of visual areas (that is, maps whose values reflect the orientation of the visuotopic gradient vector relative to the vector that points ‘geodesically’ towards V1, see Supplementary Methods 6.5).

The classifier analyses were conducted using a standard machine learning train/test/validation approach. The classifier was trained using the 210P subjects and tested against overfitting using the extra 29T subjects. The 210V subjects were used as the validation sample, and thus were not involved in the classifier training, testing, or the parcellation itself, and also shared no family relationships with the 210P or the 29T groups. A short initial run of the classifier was used to identify features that the classifier was particularly sensitive to for each area (see below and Supplementary Methods 6.6). These features were compared in each individual subject with the group average pattern to exclude subjects that were potentially misaligned with the typical subject in this region (and hence for which the group defined training labels were likely inaccurate). This area-specific set of subjects in the 210P and 29T groups were excluded from the final classifier training of each area. The classifier’s output (ranging from 0 to 1) represents the likelihood that a given vertex in a subject is part of the area being classified or part of the surrounding cortex of the searchlight. Once the classifier training weights have been generated, it is possible to classify any subject who has the 112 multi-modal maps computed, including those whose areas are misaligned with the group (see Supplementary Results and Discussion 1.4).

The trained classifier was applied to the 449 subjects and 27 repeat subjects to generate individual subject likelihood maps for each of the 180 areas in each hemisphere. These probability maps were combined by finding the largest probability for each vertex and then regularized within local neighbourhoods (see Supplementary Methods 6.7) to make an individual subject ‘winner-take-all’ parcellation. An area was considered to have been detected in a subject for the purposes of the areal detection measures (the overall classifier areal detection rate and the maps of areal detection rate for each area) if its size was between 1/3 and 3 times the size of the original population-based parcel (a pragmatic threshold chosen prior to performing the analysis that tolerates modestly greater neuroanatomical variability across subjects than the empirical range reported in cytoarchitectonic studies38,39). Probabilistic maps of each area were then created by separately averaging the individual subject winner-take-all parcellation areas for the 210P and 210V subject groups. A group maximum probability map (MPM) parcellation was then created by assigning the identity of the maximum areal probability to each vertex. The reproducibility of the parcellation was assessed by correlating these two MPM maps and by computing a Dice coefficient. In both cases the parcellation was first turned into 180 concatenated binary ROIs per hemisphere (each area was represented by a separate map, ~30,000 vertices per hemisphere, with ones for all vertices inside the area and zeros for all vertices outside). The reproducibility of the individual subject hard parcellation maps was assessed similarly. For more details, see Supplementary Methods 7.2.

Multi-modal areal fingerprints learned by the classifier were visualized using a classifier sensitivity metric. This metric was the partial derivative with respect to each feature of each area multiplied by the gradient magnitude of the feature (see Supplementary Methods 6.8). The measure indicates which areal features the classifier finds most informative when classifying a given area and whether increases or decreases in the value of the feature make the area more likely to be present. The sensitivity metric can be visualized both at the dense (vertex-wise) level for each feature and each area, or summarized at a parcel level. For each feature, the sensitivity metric was summarized at the parcel level by taking the maximum absolute value of the metric (finding the border where the feature was most influential) and using this maximum to represent the area in a parcellated or a matrix view, as shown in Supplementary Fig. 12 of Supplementary Results and Discussion 1.6.

Data reporting

No statistical methods were used to predetermine sample size. The experiments were not randomized. The investigators were not blinded to allocation during experiments and outcome assessment.

Supplementary Material

Supplementary information is available in the online version of the paper.

Extended Supplementary Methods

Neuroanatomical Supplementary Results

Supplementary Information Guide

Supplementary Results and Discussion

Acknowledgements

We thank the members of the WU-Minn-Ox HCP Consortium for invaluable contributions to data acquisition, analysis, and sharing and E. Reid and S. Danker for assistance with preparing the manuscript. Supported by NIH F30 MH097312 (M.F.G.), ROIMH-60974 (D.C.V.E.), NIH F30 MH099877 (C.D.H.), the Human Connectome Project grant (1U54MH091657) from the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research, and the Wellcome Trust Strategic Award 098369/Z/12/Z (S.M.S., J.A., C.F.B., M.J.).

Footnotes

Contributed by

Author contributions M.F.G. and D.C.V.E. designed the study and carried out the analyses. M.F.G., T.S.C., E.C.R., C.D.H., J.H., E.Y., K.U., J.A., C.F.B., M.J., and S.M.S. contributed novel methods. M.F.G., T.S.C., E.C.R., C.D.H., E.Y., J.A., C.F.B., M.J., S.M.S., and D.C.V.E. wrote the paper.

Author Information Reprints and permissions information is available at www.nature.com/reprints. The authors declare no competing financial interests. Readers are welcome to comment on the online version of the paper.

Reviewer Information Nature thanks R. Poldrack, F. Tong and T. Yeo for their contribution to the peer review of this work.

References

1. Brodmann K. Vergleichende Lokalisationslehre der Grosshirnrinde in ihren Prinzipien dargestellt auf Grund des Zellenbaues. J.A. Barth; 1909. [Google Scholar]Garey LJ, translator. Brodmann’s Localization in the Cerebral Cortex. Smith Gordon; 1994. [Google Scholar]

2. Felleman DJ, Van Essen DC. Distributed hierarchical processing in the primate cerebral cortex. Cereb Cortex. 1991;1:1–47. [PubMed] [Google Scholar]

3. Nieuwenhuys R. The myeloarchitectonic studies on the human cerebral cortex of the Vogt–Vogt school, and their significance for the interpretation of functional neuroimaging data. Brain Struct Funct. 2013;218:303–352. [PubMed] [Google Scholar]

4. Van Essen DC, Glasser MF, Dierker DL, Harwell J, Coalson T. Parcellations and hemispheric asymmetries of human cerebral cortex analyzed on surface-based atlases. Cereb Cortex. 2012;22:2241–2262. [PMC free article] [PubMed] [Google Scholar]

5. Glasser MF, et al. The minimal preprocessing pipelines for the Human Connectome Project. Neuroimage. 2013;80:105–124. [PMC free article] [PubMed] [Google Scholar]

6. Smith SM, et al. Resting-state fMRI in the Human Connectome Project. Neuroimage. 2013;80:144–168. [PMC free article] [PubMed] [Google Scholar]

7. Uğurbil K, et al. Pushing spatial and temporal resolution for functional and diffusion MRI in the Human Connectome Project. Neuroimage. 2013;80:80–104. [PMC free article] [PubMed] [Google Scholar]

8. Van Essen DC, et al. The WU-Minn Human Connectome Project: an overview. Neuroimage. 2013;80:62–79. [PMC free article] [PubMed] [Google Scholar]

9. Glasser MF, Goyal MS, Preuss TM, Raichle ME, Van Essen DC. Trends and properties of human cerebral cortex: correlations with cortical myelin content. Neuroimage. 2014;93:165–175. [PMC free article] [PubMed] [Google Scholar]

10. Glasser MF, Van Essen DC. Mapping human cortical areas in vivo based on myelin content as revealed by T1- and T2-weighted MRI. J Neurosci. 2011;31:11597–11616. [PMC free article] [PubMed] [Google Scholar]

11. Barch DM, et al. Function in the human connectome: task-fMRI and individual differences in behavior. Neuroimage. 2013;80:169–189. [PMC free article] [PubMed] [Google Scholar]

12. Caspers S, Eickhoff SB, Zilles K, Amunts K. Microstructural grey matter parcellation and its relevance for connectome analyses. Neuroimage. 2013;80:18–26. [PMC free article] [PubMed] [Google Scholar]

13. Schleicher A, Amunts K, Geyer S, Morosan P, Zilles K. Observer-independent method for microstructural parcellation of cerebral cortex: a quantitative approach to cytoarchitectonics. Neuroimage. 1999;9:165–177. [PubMed] [Google Scholar]

14. Robinson EC, et al. MSM: a new flexible framework for multimodal surface matching. Neuroimage. 2014;100:414–426. [PMC free article] [PubMed] [Google Scholar]

15. Zilles K, Amunts K. Centenary of Brodmann’s map—conception and fate. Nat Rev Neurosci. 2010;11:139–145. [PubMed] [Google Scholar]

16. Cohen AL, et al. Defining functional areas in individual human brains using resting functional connectivity MRI. Neuroimage. 2008;41:45–57. [PMC free article] [PubMed] [Google Scholar]

17. Kolster H, Peeters R, Orban GA. The retinotopic organization of the human middle temporal area MT/V5 and its cortical neighbors. J Neurosci. 2010;30:9801–9820. [PMC free article] [PubMed] [Google Scholar]

18. Wang L, Mruczek RE, Arcaro MJ, Kastner S. Probabilistic maps of visual topography in human cortex. Cereb Cortex. 2015;25:3911–3931. [PMC free article] [PubMed] [Google Scholar]

19. Gordon EM, et al. Generation and evaluation of a cortical area parcellation from resting-state correlations. Cereb Cortex. 2016;26:288–303. [PMC free article] [PubMed] [Google Scholar]

20. Shen X, Tokoglu F, Papademetris X, Constable RT. Groupwise whole-brain parcellation from resting-state fMRI data for network node identification. Neuroimage. 2013;82:403–415. [PMC free article] [PubMed] [Google Scholar]

21. Yeo BT, et al. The organization of the human cerebral cortex estimated by intrinsic functional connectivity. J Neurophysiol. 2011;106:1125–1165. [PMC free article] [PubMed] [Google Scholar]

22. Hopf A. Uber die Verteilung myeloarchitektonischer Merkmale in der Stirnhirnrinde beim Menschen. J Hirnforsch. 1956;2:311–333. [PubMed] [Google Scholar]

23. Van Essen DC, Glasser MF. In vivo architectonics: a cortico-centric perspective. Neuroimage. 2014;93:157–164. [PMC free article] [PubMed] [Google Scholar]

24. Olman CA, et al. Layer-specific fMRI reflects different neuronal computations at different depths in human V1. PLoS One. 2012;7:e32536. [PMC free article] [PubMed] [Google Scholar]

25. Polimeni JR, Fischl B, Greve DN, Wald LL. Laminar analysis of 7T BOLD using an imposed spatial activation pattern in human V1. Neuroimage. 2010;52:1334–1346. [PMC free article] [PubMed] [Google Scholar]

26. Yacoub E, Harel N, Ugurbil K. High-field fMRI unveils orientation columns in humans. Proc Natl Acad Sci USA. 2008;105:10607–10612. [PMC free article] [PubMed] [Google Scholar]

27. Zimmermann J, et al. Mapping the organization of axis of motion selective features in human area MT using high-field fMRI. PLoS One. 2011;6:e28716. [PMC free article] [PubMed] [Google Scholar]

28. Smith SM, et al. Functional connectomics from resting-state fMRI. Trends Cogn Sci. 2013;17:666–682. [PMC free article] [PubMed] [Google Scholar]

29. Hacker CD, et al. Resting state network estimation in individual subjects. Neuroimage. 2013;82:616–633. [PMC free article] [PubMed] [Google Scholar]

30. Tavor I, et al. Task-free MRI predicts individual differences in brain activity during task performance. Science. 2016;352:216–220. [PMC free article] [PubMed] [Google Scholar]

31. Van Essen DC, et al. The brain analysis library of spatial maps and atlases (BALSA) database. Neuroimage. 2016 doi: 10.1016/j.neuroimage.2016.04.002. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

32. Hill J, et al. A surface-based analysis of hemispheric asymmetries and folding of cerebral cortex in term-born human infants. J Neurosci. 2010;30:2268–2276. [PMC free article] [PubMed] [Google Scholar]

33. Van Essen DC, Dierker DL. Surface-based and probabilistic atlases of primate cerebral cortex. Neuron. 2007;56:209–225. [PubMed] [Google Scholar]

34. Filippini N, et al. Distinct patterns of brain activity in young carriers of the _APOE-ε_4 allele. Proc Natl Acad Sci USA. 2009;106:7209–7214. [PMC free article] [PubMed] [Google Scholar]

35. Abdollahi RO, et al. Correspondences between retinotopic areas and myelin maps in human visual cortex. Neuroimage. 2014;99:509–524. [PMC free article] [PubMed] [Google Scholar]

36. Griffanti L, et al. ICA-based artefact removal and accelerated fMRI acquisition for improved resting state network imaging. Neuroimage. 2014;95:232–247. [PMC free article] [PubMed] [Google Scholar]

37. Salimi-Khorshidi G, et al. Automatic denoising of functional MRI data: combining independent component analysis and hierarchical fusion of classifiers. Neuroimage. 2014;90:449–468. [PMC free article] [PubMed] [Google Scholar]

38. Caspers S, et al. The human inferior parietal lobule in stereotaxic space. Brain Struct Funct. 2008;212:481–495. [PubMed] [Google Scholar]

39. Malikovic A, et al. Cytoarchitectonic analysis of the human extrastriate cortex in the region of V5/MT+: a probabilistic, stereotaxic map of area hOc5. Cereb Cortex. 2007;17:562–574. [PubMed] [Google Scholar]