Parcellations and Hemispheric Asymmetries of Human Cerebral Cortex Analyzed on Surface-Based Atlases

Abstract

We report on surface-based analyses that enhance our understanding of human cortical organization, including its convolutions and its parcellation into many distinct areas. The surface area of human neocortex averages 973 cm² per hemisphere, based on cortical midthickness surfaces of 2 cohorts of subjects. We implemented a method to register individual subjects to a hybrid version of the FreeSurfer “fsaverage” atlas whose left and right hemispheres are in precise geographic correspondence. Cortical folding patterns in the resultant population-average “fs_LR” midthickness surfaces are remarkably similar in the left and right hemispheres, even in regions showing significant asymmetry in 3D position. Both hemispheres are equal in average surface area, but hotspots of surface area asymmetry are present in the Sylvian Fissure and elsewhere, together with a broad pattern of asymmetries that are significant though small in magnitude. Multiple cortical parcellation schemes registered to the human atlas provide valuable reference data sets for comparisons with other studies. Identified cortical areas vary in size by more than 2 orders of magnitude. The total number of human neocortical areas is estimated to be ∼150 to 200 areas per hemisphere, which is modestly larger than a recent estimate for the macaque.

Introduction

Human cerebral cortex, like that of other primates, contains a complex mosaic of areas that differ in structure, function, and connectivity. Based on postmortem architectonic analyses, many competing parcellation schemes of human cortex have been described (see Zilles and Amunts 2010). These architectonic approaches are complemented by a growing arsenal of noninvasive neuroimaging methods that use magnetic resonance imaging (MRI) to map the arrangement of cortical areas. Functional MRI (fMRI) provides information about cortical parcellation based on retinotopic maps, task-evoked activation patterns, and gradients in functional connectivity revealed by resting-state fMRI (Grill-Spector and Malach 2004; Cohen et al. 2008; Nelson et al. 2010; Wandell and Winawer 2010). Diffusion imaging and tractography enable cortical parcellation based on differences in structural connectivity (Johansen-Berg and Rushworth 2009; Mars et al. 2011). However, despite much progress, an accurate consensus cortical parcellation is lacking for the majority of human neocortex. A major reason is that differences between neighboring areas are often subtle. Major additional confounds arise from the complexity and variability of human cortical convolutions and from individual variability in the size, shape, and location of each cortical area.

Surface-based visualization and analysis offer many advantages for dealing with the complexities of human cortical convolutions (Van Essen et al. 1998; Fischl, Sereno, and Dale 1999). However, comparison of results obtained in different individuals is hampered by the dramatic variability in the pattern of human cortical convolutions and in the location of cortical areas relative to these folds. Surface-based registration (SBR) offers an attractive general approach for addressing these problems by aligning individuals to an atlas target (Van Essen 2005; Van Essen and Dierker 2007). However, SBR methods currently in widespread use rely only on local shape features and do not compensate for variability in the location of areas relative to cortical folds (Fischl, Sereno, Tootell, et al. 1999; Van Essen 2005; Goebel et al. 2006; Joshi et al. 2007; Lyttelton et al. 2007; Yeo et al. 2010; see Discussion).

The first surface-based human atlases were reconstructed from individual brains (Van Essen and Drury 1997; Van Essen et al. 1998; Van Essen 2002). However, the idiosyncrasies of cortical folding cause undesirable biases when using any particular individual as a target atlas. Such biases can be avoided by using surface-based atlases based on a population average, analogous to the population-average volumetric atlases that are widely used in human neuroimaging. The present study builds upon 2 existing surface-based atlases (FreeSurfer's “fsaverage” and Caret's “PALS-B12” atlases) and uses 2 SBR approaches: Energy-SBR and Landmark-SBR. The fsaverage atlas was generated independently for the left and right hemispheres of 40 adults using an energy minimization algorithm (Energy-SBR) that aligns continuously variable shape features (“average convexity”) in the individual and atlas surfaces (Fischl, Sereno, Tootell, et al. 1999; Desikan et al. 2006). The PALS-B12 atlas (Van Essen 2005) was generated from 12 adults using Landmark-SBR to align explicit landmark contours on the individual and atlas surfaces and to bring the left and right hemispheres into geographic correspondence with one another. Here, we registered the left and right fsaverage atlas surfaces to a common target (the “fs_LR” hybrid atlas) using a recently implemented Landmark-SBR algorithm (Van Essen et al. 2011). We also registered the fs_LR and the PALS-B12 atlases to one another, thereby enabling bidirectional migration of data between atlases.

We used these atlas refinements to address 3 general issues. 1) An analysis of shape features in the left and right hemispheres demonstrates important symmetries in population-average cortical folding patterns, quantifies the 3D positional asymmetries between the 2 hemispheres, and reveals an intriguing pattern of local asymmetries in surface area. 2) An analysis of different SBR methods confirms that the pattern of local distortions associated with SBR depends on both the algorithm used to mediate the transformations and the shapes of the individual and atlas surfaces. It demonstrates limitations of current SBR methods in compensating for individual variability. This extends and complements the comparison of a 26-landmark version of Landmark-SBR to the FreeSurfer and Brain Voyager versions of Energy-SBR carried out by Pantazis et al. (2010). 3) Numerous cortical parcellations have been mapped onto human atlas surfaces and made freely available. This enables comparisons to be made on a common atlas framework and facilitates identification of discrepancies as well as commonalities among different parcellation schemes. Informative examples illustrated herein include areas of orbital and medial prefrontal cortex (OMPFC); the middle temporal (MT) area and nearby areas in the superior temporal sulcus; and areas in medial and posterior parietal cortex.

Another objective of this study was to estimate the total number of cortical areas in humans. This number represents a fundamental aspect of brain organization in any mammalian species, yet has proven very difficult to assess even in intensively studied species like the macaque. Published parcellations that span most or all of the cortical sheet span a wide range in the total number of areas. Here, we generated a composite parcellation in which the areas identified in each region are derived from a parcellation considered to be especially reliable in that region. We estimate that human cortex contains ∼150 to 200 cortical areas. We also consider whether the number might be smaller if the heterogeneity of some regions involves fluctuations and gradients in functional and connectional properties that do not reflect well-defined cortical areas (see Discussion).

Materials and Methods

Surface-based analysis and visualization were carried out using Caret v5.616 and earlier versions (Van Essen et al. 2001; Caret, http://brainvis.wustl.edu/wiki/index.php/Caret:About). Some data sets were processed using FreeSurfer (http://surfer.nmr.mgh.harvard.edu/fswiki). This section describes the sources of data plus the methods for 1) mapping individual hemispheres to the PALS-B12 and fsaverage atlases, 2) generating the fs_LR atlases, 3) interatlas registration between PALS-B12 and fs_LR atlases, 4) surface area analyses, and 5) processing of probabilistic architectonic maps.

Data Sources

Analyses were carried out on existing data sets obtained from multiple sources. The fsaverage atlas surfaces were downloaded from http://surfer.nmr.mgh.harvard.edu/fswikialong with probabilistic architectonic maps of 13 cortical areas. The fsaverage atlas is based on 40 adults (10 each of young, middle-aged, and elderly adults, plus 10 Alzheimer's disease patients (Desikan et al. 2006).

Two independent data sets were used for analyses of cortical shape and surface area. The “OASIS-24” data set of 24 young adults (12 females and 12 males) downloaded from the Open Access Series of Imaging Studies (OASIS) archive (Marcus, Wang, et al. 2007; http://www.oasis-brains.org/) included FreeSurfer-generated surfaces and thickness maps in their native (preregistration) dimensions generated using FreeSurfer (v. 4.5) in an Extensible Neuroimaging Archive Toolkit (XNAT) processing pipeline (Marcus, Olsen, et al. 2007). The “Conte-69” data set was based on surfaces generated from 69 individual subjects (38 males and 31 females; age range 9–45) scanned under the auspices of the Conte Center for the Neuroscience of Mental Disorders. This data set has also been used in an analysis of myelin maps (Glasser and Van Essen 2011; Glasser et al. 2011).

Additional data sets were used to map cortical parcellations to atlas surfaces. These include FreeSurfer-generated surfaces and parcellations from the studies of Kolster et al. (2010), Swisher et al. (2007), Pitzalis et al. (2006), and Tsao et al. (2008); surfaces from Brewer et al. (2005) derived using mrGray (http://white.stanford.edu/∼brian/mri/segmentUnfold.htm); and Caret-generated prefrontal surface reconstructions and parcellations from Ongur et al. (2003). Additional parcellations from Brodmann (1909), Hadjikhani et al. (1998), and Press et al. (2001) had previously been mapped to the PALS-B12 atlas (Van Essen 2004; Van Essen et al. 2005). Volumetric population-average architectonic maps were obtained from the Statistical Parametric Mapping (SPM) Anatomy Toolbox (Eickhoff et al. 2005; http://www2.fz-juelich.de/inm/inm-1/spm_anatomy_toolbox/).

Registering Individual Hemispheres to Atlas Surfaces

Registration to the PALS-B12 atlas requires generation of the “Core 6” landmarks on each individual hemisphere. Landmarks were either drawn manually (http://brainvis.wustl.edu/help/landmarks_core6/landmarks_core6.html/) or by an “Automated Landmark Identification” process that includes quality control steps to review and edit landmarks (Anticevic et al. 2011; http://brainvis.wustl.edu/wiki/index.php/Caret:Operations/SurfaceBasedRegistration#Automatic_Landmarks). Once landmarks were drawn, registration to the PALS-B12 atlas landmarks was done using the Landmark Pin and Relax algorithm (Van Essen 2005).

The partial hemisphere prefrontal cortex surface reconstructions from Ongur et al. (2003) were registered to the PALS-B12 atlas surface by a process similar to that used for macaque partial hemisphere surfaces (Fig. 5 in Van Essen et al. 2011). The main difference is that landmarks were projected directly from a flat map to a partial spherical surface rather than from an inflated partial hemisphere. Registration to the PALS-B12 atlas in these cases was carried out using the Landmark Vector Difference (LVD) algorithm (Van Essen et al. 2011).

Registration of individual FreeSurfer-generated surfaces to the fsaverage atlas was carried out as part of the “recon-all” process in FreeSurfer. Figure 1 illustrates key stages in the registration of an individual right hemisphere to the fsaverage atlas and shows the spatial distortions that occur during the process. Figure 1A shows the FreeSurfer-generated white and pial surfaces plus the Caret-generated midthickness surface. A spherical map for the individual subject (Fig. 1B) includes shape features of the individual indicated by a map of average convexity (“rh.sulc,” which is analogous to the Caret-derived sulcal depth). The areal distortion map between the sphere and the midthickness surface (Fig. 1C) includes patches of local compression or expansion in regions containing high intrinsic curvature. FreeSurfer's Energy-SBR aligns the individual's average convexity map (Fig. 1B) to the population-average convexity map (“rh.avg_sulc”) (Fig. 1D), resulting in a deformed average convexity map (Fig. 1E). Comparison of Figure 1D,E demonstrates that shape features in the deformed individual map are generally well matched to those of the population-average map. However, this entails large local distortions in regions where individual and population-average shape features differ markedly (Fig. 1B vs. Fig. 1D). These registration-related distortions (Fig. 1F, ratio of spherical postregistration vs. preregistration surface tiles) approach 3-fold in many regions (yellow and green). In general, Energy-SBR distortion maps show greater local variations than those generated using Landmark-SBR (cf. Fig. 4D in Van Essen, Glasser, et al. 2011; Fig. 5D in Van Essen 2005).

A comparison of geographic correspondences between the individual midthickness surface resampled to the fsaverage mesh (Fig. 1G) and the fsaverage average midthickness surface (Fig. 1H) reveals important aspects of the FreeSurfer registration process. Highlighted nodes in the individual generally match to corresponding gyral/sulcal features in the fsaverage surface. However, in some places, this distorts local 3D relationships. For example, nodes 5 and 6 are separated by a large sulcal crease in the individual but by a gyral ridge in the fsaverage surface. Large local distortions are most likely to occur in regions of high individual variability but are not otherwise consistent across hemispheres (see Results). A map of areal distortion between the individual and the fsaverage midthickness surface (Fig. 1I displayed on the very inflated atlas surface) shows hotspots distributed across many regions, especially in regions of high individual variability in lateral prefrontal, parietal, and temporal cortex and including the region separating highlighted nodes 5 and 6.

Registration to the fs_LR Atlas

Figure 2A shows the fsaverage left and right hemisphere midthickness surfaces viewed concurrently in Caret. Although the 2 surfaces contain the same number of nodes, they are not in geographic correspondence because any given node in the left hemisphere (the “fs_L” mesh, see Table 1) corresponds to a right hemisphere node (the “fs_R” mesh) at a different location, often very distant. For example, node 1 (green) lies along the midline in occipital cortex of the right hemisphere but is much more lateral in the left hemisphere; node 2 (blue) lies on the precentral gyrus of the left hemisphere but on the superior frontal sulcus in the right hemisphere.

To bring the left and right fsaverage surfaces into geographic correspondence, both hemispheres were registered to an unbiased left–right hybrid target (Fig. 2D) using Landmark-SBR. A total of 55 landmark contours were delineated, representing corresponding features identifiable in each average midthickness surface (Fig. 2B). In general, there was little ambiguity regarding what constitutes corresponding geographic features in the 2 hemispheres. There was typically some uncertainty in specifying landmark terminations, but this was generally a small fraction of overall landmark contour length. Landmarks were projected to the corresponding left and right spherical surfaces (Fig. 2C). An average of the right plus mirror-flipped left spherical landmark contours provided an unbiased target for Landmark-SBR using the LVD algorithm (Van Essen, Glasser, et al. 2011). The spherical mesh used as the substrate for this contained the same number of nodes as the left and right fsaverage surfaces but has a more uniform node spacing in order to facilitate analyses of cortical connectivity that will benefit from the regularity of node spacing. It was created using an algorithm (“Polyhedron New”) that starts with an icosahedral mesh and subdivides it into a mesh with nodes that are nearly equally spaced and has tiles nearly equal in surface area.

After registration to the hybrid mesh, the fsaverage left and right surfaces are in good geographic correspondence, as shown by the 5 highlighted nodes on the midthickness and inflated surfaces (Fig. 2E). For example, node “2” lies on the supramarginal gyrus (SMG) in both hemispheres, even though it is in a region of pronounced hemispheric asymmetry (see Results). To reflect their distinctness despite the same node number, the underlying surface meshes are designated as “164k_fs_L” and “164k_fs_R” for the left and right hemispheres, respectively (Table 1; fs_L and fs_R for short). The hybrid atlas is designated as the “164k_fs_LR” mesh (fs_LR for short).

Deformation maps and associated data files (cf. Van Essen, Glasser, et al. 2011) that encode the fs_L-to-fs_LR and fs_R-to-fs_LR transformation are available and can be used to register FreeSurfer data sets to the fs_LR atlas. For data already represented on the fs_L or fs_R mesh (e.g., “rh.sulc.reg” generated by FreeSurfer), this entails applying the appropriate deformation map (fs_L-to-fs_LR or fs_R-to-fs_LR) to the data set(s) of interest. For data represented on the native mesh of the individual hemisphere, the 2 successive transformations (native → fs_L or fs_R → fs_LR) can be concatenated into a single individual-to-fs_LR deformation map by a 2-step process (http://brainvis.wustl.edu/wiki/index.php/Caret:Operations/FreeSurferAndFsAverage). 1) For a right hemisphere exemplar, each node of the individual subject's native-mesh sphere.reg surface (sphere “A”) is projected to the fs_R sphere by determining its barycentric coordinates within the nearest tile in the fs_R target sphere, then “unprojected” to form sphere “B” by determining its location in the corresponding tile of the fs_R-to-fs_LR deformed sphere. 2) An individual-to-fs_LR deformation map is generated by determining the barycentric coordinates between each sphere B node and the nearest tile in the target fs_LR sphere.

Interatlas Registration (fs_LR to PALS-B12)

Figure 3A shows the PALS-B12 left and right average midthickness surfaces (Van Essen 2005) displayed on the “74k_pals” mesh (Table 2). The similarities in location of the 5 pairs of highlighted nodes indicate good geographic correspondence except for modest misalignment in lateral temporal cortex (nodes 1 and 5, for details, see Figure Legend). The PALS-B12 surfaces are more irregular, especially in lateral parietal–occipital cortex. This is because only 12 subjects were averaged and also because the Landmark-SBR algorithm does not force sharp alignment of local folds in regions distant from the landmarks.

The 55 landmarks used to register the fsaverage left and right hemispheres are also recognizable on the PALS-B12 average midthickness surfaces (Fig. 3B). These landmarks were used for bidirectional registration between the fs_LR sphere (Fig. 3C same as Fig. 2C) with the PALS-B12 left hemisphere sphere (Fig. 3B) and separately with the PALS-B12 right hemisphere. However, this does not preclude possible systematic biases in the mapping between atlases, owing to differences in how shape features in the population-average atlas surfaces are related to the statistical distribution of sulci and gyri in individual subjects (see Fig. 9).

The transition between cortex and the noncortical “medial wall” region on the PALS-B12 atlas was revised (relative to Van Essen 2005) so that cortex includes all of neocortex plus pyriform cortex and transitional cortex in the hippocampal sulcus (subiculum) and basal forebrain but excluded hippocampus and other allocortical regions.

The quality of the resultant interatlas registration is high. For example, Figure 3D shows 5 highlighted nodes in good geographic correspondence on all 4 atlas surfaces. Each highlighted node is in a region of high folding variability and nodes 2 and 3 are in a region of known hemispheric asymmetry in parietotemporal cortex.

Deformation map files (cf. Van Essen, Glasser, et al. 2011) that encode the pals_L-to-fs_LR and pals_R-to-fs_LR transformations are available and can be used for mapping existing data on the PALS-B12 atlas to fs_LR. Data sets that are represented on the native mesh of a hemisphere that has already been registered to PALS-B12 can be mapped to the fs_LR atlas in a single step by a concatenation process similar to that described above for individual FreeSurfer hemisphere data sets.

Surface Area Analyses

For the OASIS-24 data set, the existing FreeSurfer-generated pial and white surfaces were converted to Caret data format and averaged to create left and right cortical midthickness surfaces for each individual. The noncortical medial wall as defined on the PALS-B12 atlas surface was mapped to each individual using a registration process that includes semiautomated landmark identification (Anticevic et al. 2011). Cortical surface area was computed for each midthickness surface after exclusion of the medial wall. To determine the impact of various atlas transformations, surface areas were recalculated for surfaces in the FreeSurfer-generated Montreal Neurological Institute (MNI) space (using talairach.xfm) and after transformation to 711-2 space by a standard MNI to 711-2 affine transformation.

For the Conte-69 data set, surface areas were computed on the cortical midthickness surfaces in native space and on the native mesh. The noncortical medial wall was excluded by discounting regions in which FreeSurfer's thickness maps had a value of zero. To create population-average surfaces for the Conte-69 data set, linear volumetric registration between the individual subject and the MNI152_T1_1mm.nii.gz was performed using the Centre for Functional Magnetic Resonance Imaging at Oxford University Linear Image Registration Tool (FLIRT). (This target is the nonlinearly derived template distributed with FSL version 4.1.7.) The resultant affine transform was applied to the FreeSurfer white, pial, and midthickness surfaces before they had been resampled to the fs_LR mesh. The individual-subject midthickness surfaces were averaged separately for the left and right hemispheres. In addition, the distortion between individual and population-average midthickness surfaces was computed (as in Fig. 1I), then averaged across all 69 individuals separately for the left and right hemispheres. The resultant average distortion maps enable surface areas computed for any region of interest on the population-average midthickness surface to be adjusted to reflect its average surface area in individual subjects, as previously demonstrated for the PALS-B12 atlas (Van Essen 2005).

To assess hemispheric asymmetries in surface area, a map of each individual's 164k_fs_LR standard-mesh surface area was computed (node i's area is the mean of its adjoining tiles' areas). For visualization, surface area maps were averaged across all subjects separately for the left and right hemispheres. The ratio of these maps yielded a population-average asymmetry map. To test for statistical significance, the individual left and right area maps were used as input to a paired t-test and 2500 iterations of surface-based Threshold-Free Cluster Enhancement (TFCE) (Hill et al. 2010; see also Smith and Nichols 2009).

Probabilistic Architectonic Maps

Surface-based probabilistic maps were available for some cortical areas, based on many individual-subject parcellations aligned to the atlas using SBR (Fischl et al. 2008) and http://surfer.nmr.mgh.harvard.edu/fswiki. These were converted to composite maps showing the most likely areal boundaries and overall extent of each area. Thresholding of the probabilistic map at a single threshold value (e.g., 0.5) was unsatisfactory for generating these composite maps because the consistency of intersubject alignment varied markedly across areas. Instead, in regions where probabilistic maps of adjacent areas overlapped, a common boundary was drawn along regions of equal probability. In regions lacking a boundary with an adjacent architectonic area, the boundary was drawn so that the size of each area on the composite map approximated the average size of the contributing individual areas.

Volumetric probabilistic maps of 44 cortical areas were obtained from the SPM Anatomy Toolbox (http://www.fz-juelich.de/inm/inm-1/spm_anatomy_toolbox; Eickhoff et al. 2005), using data sets generated by observer-independent architectonic parcellation (see above) from the Zilles and Amunts lab. These data sets had been registered from individual postmortem MR scans to the “Colin27” individual brain in MNI152 space by nonlinear volume registration (Eickhoff et al. 2005; Schleicher et al. 2005; 2009). The Colin27 MR volume was registered to MNI152_T1_1mm.nii.gz by a linear registration stage using FLIRT followed by a nonlinear stage using FMRIB's Non-linear Image Registration Tool (FNIRT). The fsaverage average MR volume was also linearly and nonlinearly registered to the MNI152 template. This transformation was inverted and concatenated to the linear transforms and the other nonlinear transform to produce a Colin27-to-fsaverage nonlinear transform. This transform was applied to the probabilistic volumetric maps of the SPM Anatomy Toolbox using nearest neighbor resampling to maintain the discrete probability values present in the original data. The center of gravity of each map was calculated and mapped to the fs_LR left and right midthickness surfaces as a stereotactic focus (coordinate).

Geodesic Distances

Geodesic distances along the cortical surface were used to compare results from different registration algorithms. These were computed by using Dijkstra's algorithm, after adding additional paths between nodes where following the surface of 2 triangles gives a shorter distance than following 2 edges. This was done in order to reduce “jaggedness” of the trajectories followed, giving a more accurate geodesic distance. In order to estimate typical intracortical geodesic distances, distances were computed along 5 representative individual midthickness surfaces rather than the fsaverage midthickness.

Data Availability

Data sets associated with this study are freely available in the SumsDB database (http://sumsdb.wustl.edu/sums/directory.do?id=8286149&dir_name=HUMAN_ATLAS_CC11). Each figure showing surface representations in this study is associated with a “scene” file that allows recreation of all of the constituent figure panels. Additional scene files provide useful starting points for various analyses. To facilitate cross-platform data migration, most surface-related files are in GIFTI format (http://www.nitrc.org/projects/gifti). All data sets can be viewed online using WebCaret without downloading data or software and can be downloaded for offline analysis using Caret or other GIFTI-compliant software.

Results

The first part of the Results includes a quantitative analysis of size and shape, including symmetries and asymmetries between the left and the right hemispheres. The second part presents a set of cortical parcellations mapped to the fs_LR atlas. We illustrate several insights gained from cross-study comparisons on a common atlas substrate and consider some of the challenges that arise when comparing data registered by different methods (Volume-Based Registration [VBR], Energy-SBR, and Landmark-SBR).

Cortical Surface Area and Variability

Table 2 shows values for average neocortical surface area and its variability in humans. For the 69 subjects in the Conte-69 cohort, average surface area is 978 ± 88 cm² per hemisphere for cortical midthickness surfaces in their native dimensions (prior to atlas registration) and excluding the noncortical medial wall region. The value is similar for the OASIS-24 data set (968 ± 70 cm²), the grand average (unweighted) for the 2 data sets is 973 cm².

Human cortical surface area is about 10% larger in males compared with females (11% for Conte-69 and 8% for OASIS-24) but nearly the same in the left versus right hemisphere (965 vs. 970 cm² for the OASIS-24 and 980 vs. 975 for the Conte-69 data set). Variability in surface area (standard deviation) is modest (6–8%) for both males and females in both data sets. In the OASIS data set, average cortical thickness is 2.44 ± 0.10 mm and is indistinguishable in left and right hemispheres and in males and females. Estimated neocortical volume (thickness times surface area) averages 236 ± 18 cm³ for each hemisphere. After registration of the OASIS-24 cohort to stereotaxic space, average cortical surface area is substantially larger (25% after registration to the MINI152 template and 13% for the 711-2B space). This reflects the generally larger dimensions of human brain atlas templates relative to actual brain size (Buckner et al. 2004; Lancaster et al. 2007; Hill et al. 2010).

Hemispheric Symmetries and Asymmetries

Figure 4A shows lateral views of the fsaverage midthickness surfaces with Cartesian axes overlaid. As noted already (see Materials and Methods, Fig. 2), the folding patterns are strikingly similar in the 2 hemispheres, down to the level of secondary and tertiary folds. This symmetry is especially noteworthy given 1) the high degree of individual variability in cortical folding, 2) the lack of strong left–right symmetry in the detailed convolutions for any individual subject (Ono et al. 1990; Van Essen 2005), and 3) substantial 3D positional offsets, especially in lateral temporal cortex (see below).

To assess the consistency of these observations, we carried out FreeSurfer segmentation and Energy-SBR registration on 69 subjects (the Conte-69 data set, see Materials and Methods). The resultant Conte-69 average midthickness surfaces (Fig. 4B) are very similar in shape to the fsaverage midthickness surfaces (Fig. 4A) in terms of the relative positions and trajectories of minor as well as major folds. Thus, FreeSurfer's Energy-SBR generates population-average midthickness surfaces whose shape is largely independent of the choice of individual contributing hemispheres.

Despite the symmetry in folding patterns, there are important hemispheric asymmetries in both the fsaverage and the Conte-69 data sets. One asymmetry involves 3D positional offsets of gyral and sulcal features on the average midthickness surfaces. Consistent with previous hemispheric asymmetry studies (Steinmetz et al. 1990; Habib et al. 1995; Rubens et al. 1976; Toga and Thompson 2003), the left SMG (blue arrows) is more than 1 cm posterior to the right hemisphere (blue arrows in Fig. 4A,B). A more subtle but significant shape asymmetry is evident near the ventral tip of the postcentral gyrus, where the left hemisphere has a gyral bulge that is less pronounced in the right hemisphere (white arrows). To better quantify the positional asymmetries, maps of 3D coordinate differences (distance between corresponding vertices in the right and mirror-flipped left hemisphere average midthickness surfaces) were computed and displayed on the inflated fsaverage surface for both the fsaverage (Fig. 4C) and the Conte-69 (Fig. 4D) data sets. The most prominent asymmetries are near the SMG and angular gyrus (AG) (bright red, orange, yellow, 8–14 mm in magnitude, with the left hemisphere posterior to corresponding right hemisphere nodes) and in inferior temporal cortex (left hemisphere dorsal to corresponding right hemisphere nodes). An asymmetry in the ventral calcarine sulcus largely reflects an occipital petalia, in which the left occipital pole is displaced toward the right (Toga and Thompson 2003).

Some of the positional asymmetries may be associated with asymmetries in local cortical surface area, as has been reported for the planum temporale (PT) even without the aid of explicit surface reconstructions (e.g., Habib et al. 1995). Using the Conte-69 data set mapped to the fs_LR atlas, we addressed this issue quantitatively and across the entire cortical sheet, with surprising results. We computed the average surface area associated with each tile in the left and right hemisphere Conte-69 data set (see Materials and Methods). Figure 5A shows a smoothed asymmetry map (left/right average surface area) displayed on the inflated left hemisphere. Asymmetries strongly favoring the left hemisphere (1.5- to 2-fold, yellow and bright red) are most prominent in the Sylvian Fissure (parietal operculum and anterior to Heschl's gyrus [HG]) and in medial temporal cortex (lingual gyrus and collateral sulcus). Asymmetries strongly favoring the right hemisphere (1.5- to 2-fold, green) include patches near the AG and also in dorsomedial prefrontal cortex. In many other regions, the differences are sufficiently consistent across subjects to yield high t-statistic values (|t| > 10 in many regions). A large expanse of cortex passes statistical significance (P < 0.01) by the TFCE method (Fig. 5B). The significant left-biased and right-biased regions each occupy about one-third of cortical surface area (35% and 29%, respectively) after compensating for distortions of surface area on the average midthickness surfaces averaged across the left and right hemispheres (Fig. 5C; see Materials and Methods).

Surface area asymmetries in the Sylvian Fissure are of particular interest. Figure 5D shows the same smoothed asymmetry map as in panel A but displayed on expanded (top) and full-hemisphere (bottom) views of the very inflated left hemisphere along with the registration landmarks used in this region. The asymmetries can be related to subtle but important asymmetries in the average midthickness surfaces, as seen on maps of average cortical folding on the very inflated (Fig. 5E) and average midthickness surfaces (Fig. 5F) for the left (top) and right (bottom) hemispheres. The hotspot of left-biased cortex in the parietal operculum is centered within area OP1 (see Fig. 14), between landmarks 1–4, just posterior to the aforementioned gyral ridge that is more prominent in the left hemisphere (blue asterisk in Fig. 5 same as white arrows in Fig. 4A,B). Another left-biased hotspot occupies the anterior bank of HG and adjoining insular cortex. The PT as defined by Habib et al. (1995) has only a modest (average 13%) bias for the left hemisphere even though it has been reported to be much larger in the left hemisphere (Toga and Thompson 2003; see Discussion). Near the ventral tip of the central sulcus (landmark #5 in yellow), a patch of right-biased cortex is centered on a sulcal crease that is deeper in the left hemisphere (black arrow, between landmarks 2, 5, and 6). It lies between areas 44 and OP4 (see Fig. 14). These hemispheric asymmetries are unlikely to reflect bias or inaccuracy in positioning of landmark contours, which run along prominent features in the average midthickness surfaces. While the magnitude of the asymmetry would be impacted if the landmarks were positioned slightly differently, we are confident that these hotspots represent biologically significant asymmetries (see Discussion).

Overview of Human Cortical Parcellations

Ten human cortical parcellations were registered to the fs_LR atlas. Table 3 summarizes general characteristics of each parcellation scheme and how it was registered to the atlas. It includes the study name and abbreviation (columns 1 and 2); the cortical lobe(s) analyzed (column 3) the number of reported areas (column 4); the parcellation modality, based on architectonics or retinotopic fMRI mapping (column 5); the surface reconstruction method (column 6); the registration method (column 7); the atlas used as the initial SBR target (column 8); and the figures where illustrations appear in the present study (column 9).

We first compare a set of probabilistic architectonic maps that were mapped to the fs_LR atlas using both VBR and Energy-SBR. This expands upon a previous report comparing Energy-SBR with VBR (Fischl et al. 2008). We then compare Landmark-SBR to Energy-SBR, followed by consideration of parcellations that were initially registered to the PALS-B12 atlas. Finally, we illustrate a composite human cortical parcellation based on 5 published studies.

Observer-Independent Probabilistic Architectonic Maps

The location, extent, and variability of 44 cortical areas have been charted in studies that apply observer-independent cytoarchitectonic analysis methods to postmortem human brains (see Schleicher et al. 2005, 2009). This entailed registering each postmortem MR volume to the “Colin” single-subject atlas brain using nonlinear VBR (nl-VBR, Geyer; Eickhoff et al. 2005), which were combined to generate a probabilistic volumetric map of each architectonic area. Thirteen architectonic areas were also mapped to individual-subject surface reconstructions and registered to the fs_L and fs_R atlases by Fischl et al. (2008) using Energy-SBR. These were registered to the fs_LR atlas using intraatlas Landmark-SBR (see Materials and Methods).

Figure 6 illustrates differences between nl-VBR and Energy-SBR, using area 2 as an exemplar. In individual subjects, area 2 is located largely or entirely on the anterior bank of the postcentral sulcus (PoCeS) (arrows in Fig. 6A). Owing to imperfect intersubject alignment with nl-VBR, probabilistic area 2 (Fig. 6B) spreads 30 mm along the anteroposterior axis at this slice level, spanning both banks of the PoCeS (red arrow) and encroaching on the posterior bank of the CeS (yellow arrow), part of the intraparietal sulcus (blue arrow), and even the anterior bank of the CeS. Additional VBR steps aligned the probabilistic area 2 volume to the fsaverage atlas volume (Fig. 6C; see Materials and Methods) and associated fsaverage surface (red contour in Fig. 6C). Volume-to-surface mapping enabled visualization of VBR-derived area 2 on the fsaverage atlas surface (Fig. 6D,E), revealing dispersion across neighboring sulci similar to that shown in the volume slice (Fig. 6C). In contrast, the Energy-SBR-derived map of probabilistic area 2 (Fischl et al. 2008) shows much less dispersion and is largely restricted to the anterior bank of the PoCeS (Fig. 6F). SBR also shows a higher peak probability (red = 1.0 = complete overlap) than does VBR (orange = 0.7). The VBR map occupies almost twice the surface area on the fsaverage midthickness atlas (44 vs. 24 cm²), and in some places (white nodes), it extends more than 2–3 cm beyond the SBR map when measured in terms of average geodesic distance along representative midthickness surfaces (see Materials and Methods).

Additional insights come from comparing VBR versus SBR results for area hOc5 (putative area MT, see below) in both left and right hemispheres (Fig. 7). 1) Lower alignment consistency. Maximal alignment consistency was only 0.5 (green) for both VBR (Fig. 7A) and SBR (Fig. 7B), confirming that alignment is poorer in regions of high folding variability (Fischl et al. 2008). 2) Smaller SBRversus VBR differential. VBR maps occupied greater surface area (10.2 vs. 7.7 cm² on the right, 7.9 vs. 7.6 cm² on the left), but the difference was less pronounced than for area 2. This suggests that the improved alignment by Energy-SBR versus nl-VBR is more modest in regions of high folding variability. 3) SBR outliers. Both left and right SBR maps include outliers (arrows), in which a map for a single case (P = 0.1) is largely or completely isolated from the rest of the distribution and in several instances is very elongated in shape. Such cases might arise if Energy-SBR forces large local distortions in order to match the individual subject to the population average (cf. Materials and Methods, Fig. 2). 4) Spatial offsets between VBR and SBR. The white highlighted node in the left hemisphere is centered on the peak of the left SBR map (Fig. 7C) but is at the posterior margin of the VBR map (Fig. 7A); the centers of the VBR versus SBR probabilistic maps are offset by ∼1 cm on the average midthickness surface. In the right hemisphere, the geographically corresponding node is in the posterior part of both probabilistic maps (Fig. 7B,D), and the offset between SBR and VBR maps is less pronounced. Thus, the spatial offsets between registration algorithms are not consistent across hemispheres or data sets.

Figure 8A shows probabilistic maps of 5 nonoverlapping areas (2, 4a, 45, hOc5, and 17) on lateral and medial views of the inflated right hemisphere surface. Area V1 is the most consistently aligned, with alignment consistency values exceeding 0.8 (yellow and red) over most of its extent. Areas hOc5 and 45 are least well aligned; for hOc5, the maximum alignment consistency is 0.5 for both the left and the right hemisphere maps.

A composite map of all 13 SBR-derived architectonic areas was generated for both hemispheres (Fig. 8B,C). Areal boundaries were determined by assigning areas to the highest probability in regions where probabilistic maps overlapped and by estimating the most likely extent of the typical individual (see Materials and Methods). Architectonic areas are generally symmetric in location and extent on the 2 hemispheres. Area 17 (hOc1) occupies the posterior calcarine sulcus up to its junction with the parietooccipital sulcus (about 1.5 cm from the corpus callosum and the subiculum). It is completely surrounded by area 18 (hOc2), but area 18 is only 4–5 mm wide in this region on the composite map. This contrasts with the macaque in terms of the anterior limit of area 17 and the relationship of area V2 to V1 (Van Essen et al. 1982, 2011; see Discussion).

Comparing Energy-SBR versus Landmark-SBR

It is useful to estimate the magnitude and spatial pattern of differences arising when data are registered by different SBR methods. A comparison between FreeSurfer Energy-SBR and Landmark-SBR is of particular interest, since both methods were used for mapping cortical parcellations to atlas surfaces (Table 3). We analyzed surfaces from 4 hemispheres registered to the fs_LR atlas, using Landmark-SBR and independently using Energy-SBR as the initial registration step. Figure 9 shows results for an exemplar hemisphere processed by registration to PALS-B12 (Fig. 9A) and separately by registration to fs_R (Fig. 9B). Both surfaces were subsequently mapped to fs_LR by interatlas registration (pals_R-to-fs_LR or fs_R-to-fs_LR). The 2 versions of the individual midthickness surfaces (Fig. 9A,B) are indistinguishable in shape, but each has different point-to-point correspondences with the fs_LR midthickness surface. For example, highlighted node 1 lies on the crown of the precentral gyrus in all 3 surfaces (Fig. 9A–C), indicating that the 2 methods are in correspondence at this location. In contrast, highlighted node 2 lies on the inferior frontal gyrus in all 3 surfaces, but it is 24 mm farther anterior in the fs_L-registered surface (Fig. 9B) than in the PALS-B12-registered surface (Fig. 9A). Figure 9D shows an overall map of 3D coordinate differences (distance between corresponding nodes in the 2 resampled individual and thickness surfaces) displayed on the very inflated surface. High coordinate difference values (yellow and red) are concentrated in regions of high interindividual folding variability in frontal, temporal, and parietal cortex.

The 3D difference maps for the left hemisphere of the same case (Fig. 9D, lower right) and for both hemispheres of the second case (Fig. 9E) show a similar general pattern in which large (>1 cm) 3D coordinate differences are mainly restricted to regions of high folding variability. However, there are numerous differences in the detailed patterns. The coordinate offset is less than 1 cm for most cortical regions in all 4 hemispheres examined (blue and purple regions) but it exceeds 2 cm in a few places, similar in magnitude to that reported by Pantazis et al. (2010, their Fig. 9) for different pairwise comparisons of SBR algorithms (see Discussion). This analysis provides objective data on where the results of Landmark-SBR versus Energy-SBR differ, but it does not reveal which method achieves greater consistency of intersubject alignment (see Discussion).

To explore whether there is a systematic bias between the 2 registration methods (to fs_LR via PALS-B12 vs. fs_L or fs_R), we computed the average coordinate differences separately along x-, y-, and z-axes for all 4 hemispheres (Fig. 9F–H). The population-average differences are less than 5 mm over most of cortex but in some regions are up to 15 mm along one or more axes. While this suggests potential systematic bias, replication using a larger number of hemispheres would be necessary to confirm this pattern. These findings highlight the need for caution when interpreting results that involve different SBR methods (see below).

Comparison between Retinotopic and Architectonic Areas

Numerous fMRI studies have used retinotopic mapping to chart visual areas in individual subjects. Previous efforts to validate such results have included comparisons between retinotopic V1 and architectonic area 17 (Hinds et al. 2009 also Bridge et al. 2005) and also volume-based comparisons with V2, V3, and V4 (Wohlschlager et al. 2005; Wilms et al. 2010). These studies have demonstrated significant overlap between retinotopic and architectonic areas, but the volumetric comparisons are substantially impacted by intersubject alignment errors (cf. Fig. 6). Here, we compared 3 retinotopic areas (V1, V2, and MT), all registered via PALS-B12 as the initial target (Table 4), with 3 architectonic areas (17, 18, and hOc5), all registered initially to fsaverage. Figure 10 displays retinotopic V1 (dark orange) and V2 (yellow) in 2 hemispheres, KPO10 case 1 (Kolster et al. 2010) and SHM07 case 1 (Swisher et al. 2007) in relation to the most likely boundaries of probabilistic areas 17 and 18 (black contours, see also Fig. 10E,F). The top row shows cortical areas on posterior views of the individual hemispheres (Fig. 10A,C) and the fsaverage midthickness surface (Fig. 10B,D). In KPO10 case 1 (Fig. 10A,B), retinotopic V1 and V2 boundaries are in reasonable agreement with average area 17 and 18 borders, though retinotopic V2 is narrower than the average width of architectonic area 18. SHM07 case 1 (Fig. 10C,D) shows a more pronounced mismatch between retinotopic and architectonic areal boundaries. In this case, retinotopic V1 barely extends to the occipital pole, and much of retinotopic V2 lies within the expected extent of architectonic area 17. This is unlikely to reflect a bias in registration of the individual to the atlas because on both the individual and the atlas surfaces, V1 and V2 occupy similar locations relative to the calcarine sulcus (Fig. 10C,D). Instead, this discrepancy is more plausibly attributed to inaccurate estimation of retinotopic fMRI boundaries in this subject.

In medial views (Fig. 10, middle row), the retinotopic maps of V1 and V2 are restricted to central field representations in the posterior half of areas 17 and 18, consistent with the fact that visual stimulation was restricted to central eccentricities (see Figure Legend). Retinotopic V1 lies mainly within the calcarine sulcus in both the individual and the fsaverage surfaces. Retinotopic V2 occupies a narrow strip dorsal and ventral to V1. Its relationship to architectonic area 18 most readily discerned on a posteromedial view of the very inflated atlas surface, showing retinotopic areas and architectonic boundaries for the KPO10 (Fig. 10G) and SHM07 (Fig. 10H) individual cases and for a composite of 4 hemispheres from KPO10 (Fig. 10I). The concordance between retinotopic and architectonic boundaries is reasonable in many regions. However, there are substantial discrepancies in the outer boundaries of area V2/18 for the entire SHM07 case and for parts of ventral area V2/18 in the KPO10 individual and composite.

Several factors may in principle contribute to such discrepancies. These include 1) registration inaccuracies and algorithm-dependent biases such as those already illustrated (Figs 1 and 9); 2) retinotopic border inaccuracies arising anywhere along the complex sequence of fMRI data acquisition and analysis steps, including registration of blood oxygen level–dependent echo-planar imaging volumes to the anatomical volumes (see also Van Essen et al. 2011); 3) function-folding variability, in which an apparent boundary mismatch between any individual subject and the population average may lie within the normal range of variability; and 4) different areal definitions, if the boundaries of an architectonic area (e.g., area 18) do not coincide with those of a retinotopic area (e.g., V2) when determined accurately in the same individual. In general, the relative contributions of each factor are likely to vary across cases and even across regions within any given case (see Discussion).

Architectonic area hOc5 overlaps extensively with a region of motion-sensitive cortex that includes area MT (Wilms et al. 2005; Malikovic et al. 2007). However, motion-selective fMRI activations extend across multiple cortical areas in lateral occipitotemporal cortex (Tootell et al. 1995; Huk et al. 2002; Kolster et al. 2010), and it has yet to be established whether retinotopic MT and hOc5 are coextensive. Kolster et al. (2010) used retinotopic fMRI mapping to identify area MT and other areas within this motion-sensitive region. Figure 11 compares probabilistic hOc5 with their retinotopic area MT, using data from 4 hemispheres in their study.

In the individual KPO10 hemispheres, retinotopic MT was mainly gyral in 2 hemispheres and mainly sulcal in 2 others (red in Fig. 11A,B; see also Fig. 10G,I in Kolster et al. 2010). Figure 11C shows probabilistic maps after registering MT and other KPO10 retinotopic areas to the fs_LR atlas. Only areas that overlap with the composite map of hOc5 (black contour) are displayed; this includes areas MT, FST, and V4t in both hemispheres plus MSTv on the right and LO2 and PITd on the left. Intersubject alignment for these retinotopic areas is comparable to that found for probabilistic architectonic hOc5 (Fig. 11D) even though the registration methods were different (via PALS-B12 and Landmark-SBR for the retinotopic areas and via fs_L/fs_R and Energy-SBR for hOc5). Figure 11E shows a composite map of the KPO10 retinotopic areas generated using all 4 cases and displayed on the left and right inflated atlas surfaces. By these estimates, retinotopic MT and V4t both overlap extensively with but are smaller than hOc5. However, the retinotopic MT maps are incomplete as they do not include the peripheral field representation (cf. Fig. 10 for V1 and V2 maps). The missing peripheral field representation for MT would likely be dorsal (upward) on the inflated surface and would not increase the overlap with hOc5. The misalignment between the centers of probabilistic retinotopic MT and probabilistic architectonic hOc5 may be attributable to systematic biases between the 2 registration methods in this region, especially along the z-axis (Fig. 9H). Alternatively, MT may correspond only to the dorsal subdivision of hOc5, as suggested by Malikovic et al. (2007). Overall, the hypothesis that hOc5 includes retinotopic MT is strengthened by this analysis. However, the apparent size discrepancy makes it plausible that architectonic hOc5 also includes additional visual areas (e.g., KPO10 areas V4t and/or MSTv).

Additional Retinotopic Areas

Numerous additional retinotopic areas have been charted in various fMRI mapping studies. Figure 12 shows maps from 4 such studies that have been mapped to the fs_LR atlas using Landmark-SBR via PALS-B12. These are displayed on medial (top row) dorsoposterior (bottom row) views of the inflated atlas surface. Swisher et al. (2007) mapped the location and topographic organization of 14 retinotopic areas. Their area SHM07 V3 lies almost entirely within architectonic area 18; adjoining areas V4 and LO1 lie just outside the most likely area 18 boundary in occipitotemporal cortex. SHM07 area LO2 lies posterior to hOc5 in lateral occipital cortex. In dorsal occipital and medial parietal cortex, they identified a string of 8 areas, including SHM07 V3B, V3A, V7 (a.k.a. IPS0), and IPS1–4 (Fig. 12A). Areas IPS0–IPS4 run mainly along gyral cortex in dorsomedial parietal cortex but encroach into the medial bank of the intraparietal sulcus.

Brewer et al. (2005) and Wandell et al. (2007) in combination mapped 9 retinotopic areas (Fig. 12B). Their areas BLW05 V1—V3, V3A, V3B, and V7 (IPS0) are in general agreement with the corresponding areas mapped by Swisher et al. (2007) in posterior parietooccipital cortex. In ventromedial occipitotemporal cortex, the Brewer et al. map includes BLW05 areas hV4, VO-1, and VO-2. SHM07 area V4 overlaps extensively with BLW05 hV4 but also partially with VO-1 and VO-2.

Kolster et al. (2010) mapped 10 visuotopic areas, including a visual field cluster that includes MT, MSTv, FST, and V4t, and nearby areas PITd, LO1, and LO2 (Fig. 12C; see also Fig. 11C). KPO10 areas LO1 and LO2 are in reasonable agreement with the corresponding areas of Swisher et al. As already noted, their V4t lies mainly within architectonic hOc5, and its relationship to areas in other retinotopic parcellation schemes is unclear.

Human area V6 (hV6) is a retinotopic area in medial parietal cortex that has been mapped using visual stimuli extending into the far periphery (Pitzalis et al. 2006). PGH06 area hv6 lies in the dorsal tip of the parietooccipital sulcus near the fundus (Fig. 12D), adjacent to areas V7 and IPS1 (cf. Fig. 12A,B). It lies several centimeters dorsal and anterior to areas 17 and 18, and in this respect differs markedly from the more proximate location of macaque area V6 (mV6) relative to V1 and V2 (Van Essen et al. 2011).

Human OMPFC

Ongur et al. (2003) used cytoarchitectonic, myeloarchitectonic, and chemoarchitectonic criteria to map the layout of 17 areas in OMPFC. These OMPFC areal maps are displayed on ventral views of an individual subject (Fig. 13A) and after registration to the fs_LR atlas (inflated surface, Fig. 13B). Figure 13C shows probabilistic maps for 3 exemplar OMPFC areas (4 hemispheres of 2 subjects). There is good overall alignment despite individual differences in absolute size and location relative to gyral/sulcal landmarks. Boundaries representing the most likely transition between areas were used to generate a composite parcellation (black outlines in Fig. 13D) on medial, ventral, and lateral inflated surface views.

Size Range and Estimated Number of Cortical Areas

Figure 14A shows a composite map of 52 cortical areas from 5 published parcellation schemes on the fsaverage midthickness surface. In regions where areas from competing schemes overlap, areas were selected to avoid substantial overlap (see Figure Legend). The inflated surfaces (Fig. 14B) also show the centers of gravity for an additional 31 probabilistic areas from observer-independent cytoarchitectonics that have been mapped volumetrically but not by SBR. While these additional areas span a substantial cortical domain, major portions of the temporal, frontal, and parietal lobes have yet to be charted by modern architectonic methods that can be accurately mapped to surface-based atlases (see Discussion). For reference, Figure 14C shows 43 areas from the classical Brodmann (1909) parcellation. These areas are generally larger than revealed by modern parcellations, and many of the borders are likely to be inaccurate. Nonetheless, Brodmann areal parcellations remain widely used in neuroimaging studies.

For surface area estimates, the area of each parcel on the fsaverage atlas was adjusted by a distortion–compensation factor based so that it represents surface area or a typical individual midthickness surface (see Materials and Methods). These 52 parcels had an average distortion-compensated area of 6.6 cm² and an aggregate surface area of 344 cm², which is 35% of total right hemisphere cortical surface area. As shown in a size frequency histogram (Fig. 14D), these areas vary more than 200-fold in surface area. Four areas exceed 16 cm² (areas V1, V2, 10p, and 6), whereas several areas are under 0.2 cm² (Iam, Iapm, 25, and 32pl). However, the great majority of areas (36 of 52) are within an 8-fold size range, between 2 and 16 cm².

If cortical areas in regions outside the identified areas in Figure 14 have the same average size as the 52 identified areas, it would suggest that human cortex has a total of about 148 cortical areas in each hemisphere. This would be an underestimate if most of the uncharted areas are comparable to the nonprimary and secondary and somatosensory areas. If these areas are excluded (plus Iam and Iapm), which may have been incompletely charted, then average cortical area size is 4.8 cm². This would equate to 131 areas in unassigned cortex and a total of 183 areas overall. Given the many uncertainties in this analysis, it seems reasonable to suppose that the total number is likely to be in the range of 150–200 cortical areas per hemisphere (see Discussion). This contrasts with only 43 areas shown on the panhemispheric classical Brodmann parcellation (Fig. 14C).

Discussion

Two broad sets of findings emerge from this study. A quantitative analysis of cortical shape revealed several features that are not evident from visual inspection of individual subjects. This includes hemispheric symmetries and asymmetries in the folding pattern of population-average midthickness surfaces plus an unexpected pattern of hemispheric asymmetries in cortical surface area. An analysis of cortical parcellations derived from many published studies capitalizes on recent progress in obtaining accurate surface maps of cortical areas over a substantial portion of the cortical mantle. However, a consensus parcellation is lacking for most of cortex, owing to inadequate data and/or conflicts among published parcellations, 4 general issues warrant further discussion: 1) the advantages and limitations of surface-based atlases and registration; 2) interatlas transformations and coordinate specification; 3) human hemispheric asymmetries and their functional significance; and 4) cortical parcellations within and across species.

Registration Algorithms and Atlases

Standard VBR and SBR algorithms, including those used in the present study, constrain registration using structural or geographic features as a surrogate for underlying functional organization. The fidelity of registration can be evaluated based on the alignment of geographic features themselves but also on cortical parcellations and other types of functional mapping.

For human cortex, assessing the fidelity of geographic alignment is confounded by the variability of cortical convolutions. Many geographic features (e.g., the anterior tip of the intraparietal sulcus) are ambiguous or ill defined in individual subjects, making it fundamentally problematic to determine whether different individuals are well aligned to an atlas in that region. Population-average surfaces provide a valuable anatomical substrate for assessing statistical regularities that are mostly obscured by individual variability. The 3D configuration of a population-average surface reflects 3 main factors: 1) the shapes of the contributing individual hemispheres; 2) shape features imposed by the target atlas; and 3) the irregular pattern of deformations generated by the registration algorithm when applied to each individual. In comparing 3 population-average midthickness surfaces (fsaverage, Conte-69, and PALS-B12), both left and right hemispheres of all 3 midthickness surfaces shared many shape characteristics. The fsaverage and Conte-69 average midthickness surfaces registered using Energy-SBR are especially similar (Fig. 4), indicating that the idiosyncratic shape features of contributing individuals have minimal impact if enough subjects are included. Another comparison of interest is the MNI population-average surface generated using the CIVET Energy-SBR algorithm and a different target atlas (Lyttelton et al. 2007; Im et al. 2008). This atlas shares many shape features with the PALS-B12 atlas (Lyttelton et al. 2008) and even more with the fsaverage atlas. Thus, most secondary and many tertiary shape features of population-average surfaces are consistent across target atlases, registration algorithms, and choice of contributing subjects.

Published assessments of alignment fidelity achieved in humans mainly involve comparisons between SBR and VBR. For example, Landmark-SBR achieves better alignment of identified sulci than does linear or low-dimensional nl-VBR (Van Essen 2005; Anticevic et al. 2008). Fischl et al. (2008) demonstrated that alignment consistency for probabilistic architectonic areas is much better for Energy-SBR than nl-VBR in regions of low folding variability but that the differences are smaller in regions of high variability. We confirmed this finding, showing that the differences can exceed 2 cm in geodesic distance along the cortical surface that the center of gravity of probabilistically mapped areas can differ significantly between nl-VBR and Energy-SBR (Figs 6 and 14).

Pantazis et al. (2010) carried out a detailed comparison of intersubject alignment consistency achieved by FreeSurfer, Brain Voyager, and their own 26-Landmark-SBR method. They demonstrated significant differences across methods, with Landmark-SBR performing best relative to a trained neuroanatomist. As expected, FreeSurfer and Brain Voyager aligned cortical folds (curvature) better on average, whereas Landmark-SBR aligned landmarks better. In a similar vein, we found that the alignment differences for Landmark-SBR versus Energy-SBR applied to the same hemispheres are usually modest (<1 cm in most locations) but can exceed 2 cm in regions of high variability (Fig. 9). The magnitude and direction of these differences depend on the idiosyncrasies of local distortions produced by different registration methods and are not readily predictable. There may also be systematic differences in some regions, as suggested by the population-average coordinate differences for 4 hemispheres (Fig. 9F–H).

For the various probabilistic maps registered to the fs_LR atlas, the residual variability present in these maps arises from multiple factors, including 1) inaccuracies or different criteria for identifying areal boundaries in individual subjects, 2) interindividual variability in areal size, which exceeds 2-fold even for well-defined cortical areas in humans as in the macaque (Andrews et al. 1997; Dougherty et al. 2003; Eickhoff et al. 2005; Caspers et al. 2008); 3) variability in the location of areal boundaries relative to the folds used as landmarks; and 4) local distortions associated with the registration method. The modest misalignment between the average location of architectonic hOc5 (registered via Energy-SBR) and retinotopic MT (registered via Landmark-SBR) may be attributable to systematic bias (cf. Fig. 9).

Given the limitations of existing SBR and VBR algorithms, improved registration methods are sorely needed. It is especially important that registration be constrained by information related to functional organization, not just shape (geographic) features. For example, Sabuncu et al. (2010) used task-evoked fMRI data (passive viewing of movie segments) as a registration constraint and demonstrated substantially improved alignment compared to anatomical (shape-based) Energy-SBR. This approach can be extended to include data from a broader range of fMRI tasks, but it is predicated on acquiring fMRI task data that yields robust activation patterns in each subject. An attractive alternative (or adjunct) will be to make use of architectonic features revealed by structural imaging. This includes “myelin maps” of cortical gray matter that reveal numerous cortical areas in individuals as well as population averages (Glasser and Van Essen 2011). A more ambitious approach will be to incorporate registration constraints based on connectivity patterns revealed by noninvasive imaging in individual subjects. This includes structural connectivity revealed by diffusion MRI and tractography (Johansen-Berg and Rushworth 2009) plus functional connectivity revealed by resting-state fMRI (Fox et al. 2009; Greicius et al. 2009; van den Heuvel and Hulshoff Pol 2010). “Connectome-based registration” offers the prospect of improving intersubject alignment based on a vast amount of connectivity data. However, it will require development of novel registration algorithms and it will benefit from higher quality human connectivity data generated by several ongoing efforts, including the Human Connectome Project (http://www.humanconnectome.org/consortia). Assuming that such efforts are successful, population-average atlases will continue to evolve. Rather than promoting convergence on any single atlas currently in existence, a high priority should be on developing objective methods for assessing how different algorithms and atlases perform in reducing intersubject variability and for facilitating interatlas registration.

Even as registration methods continue to improve, there will always be great value in detailed analyses of results obtained in individual subjects, along with assessment of population-average results on an atlas. One efficient strategy for capitalizing on both types of analysis is to view individual-subject data on resampled, standard-mesh versions of the appropriate individual subject's hemisphere (e.g., midthickness, pial, white, or inflated), as was done for retinotopic and architectonic maps in the present study (Figs 10–13). This allows concurrent comparisons of individual specific and atlas-related data sets in a common coordinate system

Interatlas Transformations and Coordinate Specifications

Given the multiplicity of surface-based atlases in widespread use, it is important to have efficient ways of transforming data between atlases, just as transformations are available for mapping between various volumetric atlas spaces (Lancaster et al. 2007). Interatlas transformations were achieved in the present study using Landmark-SBR and a large number (55) of identified geographic landmarks, thereby achieving high fidelity alignment. Each interatlas registration is encoded by a deformation map file that can be used to quickly map any data set of interest from one atlas to the other. Moreover, the transformations can be carried out using command-line operations, with inputs and outputs in standard GIFTI format. Extending this to other surface-based atlases such as the MNI surface-based atlas (Lyttelton et al. 2007, 2008) will further facilitate migration of data across atlases and improve cross-study comparisons.

Spatial location on a surface-based atlas can be uniquely specified in 3 distinct ways: by (x, y, z) coordinates in stereotaxic space; by latitude and longitude (θ, φ) in spherical standard coordinate space (Fischl, Sereno, and Dale 1999; Van Essen 2002, 2005); and by the node number in a standard-mesh surface. Each of these representations serves useful and complementary purposes. An advantage of including the node number is that it can immediately be linked (through machine-readable lookup tables) to the corresponding node in any other atlas or individual surface registered to the atlas. Inclusion of all 6 numbers rather than just the conventional triplet of stereotaxic coordinates in tables that report activation foci and the centers of ROIs would facilitate future cross-study comparisons.

Cortical Size, Symmetry, and Asymmetry

Our value of 1946 cm² for the average cortical surface area of the 2 hemispheres combined exceeds the values reported in some studies—1569 cm² in Van Essen and Drury (1997), 1501 cm² in Henery and Mayhew (1989), 1468–1670 cm² in Blinkov and Glezer (1968), 1678 cm² in Jouandet et al. (1989), and 1752 cm² in Klein et al. (1998). It is smaller than the value of 3031 cm² reported by Elias and Schwartz (1971) and similar to the value of 1906 cm² reported by Tramo et al. (1995). Estimates based on surfaces registered to a standard atlas are generally inflated owing to population-average atlas volumes being larger than individuals (Buckner et al. 2004).

Structural asymmetries of human cortex have been studied for more than a century. The most prominent are the asymmetric trajectories of the posterior Sylvian Fissure and the superior temporal sulcus (Eberstaller 1884; Cunningham 1892; Loftus et al. 1993; Toga and Thompson 2003; Ochiai et al. 2004; Van Essen 2005). Here, we characterized these and other positional asymmetries quantitatively by measuring the distance between corresponding points on population-average right and mirror-flipped left hemisphere surfaces. The asymmetries extend over a large swath of lateral temporal and parietal cortex, centered on the angular and supramarginal gyri. They are similar in spatial pattern and magnitude for the 2 independent populations that were analyzed (fsaverage and Conte-69).

Our analysis of asymmetries includes the first quantification of hemispheric differences in surface area across all of neocortex analyzed in a large cohort of subjects. The largest areal asymmetry is in the posterior part of area OP1 in the parietal operculum (Fig. 14). Eickhoff et al. (2006) reported the volume of OP1 to be slightly but not significantly larger in the left versus right hemisphere (6.9 vs. 5.2 cm³). Evidence for a functional asymmetry in this region comes from an fMRI study showing stronger auditory responses on the left to tool sounds versus animal vocalizations (Lewis et al. 2005; Lewis 2006) but a reversed pattern in left handers (Lewis et al. 2006). Tractography analyses using diffusion imaging data suggest an asymmetry of parietotemporal connections in this region, including connections from the supramarginal and angular gyri to the middle temporal gyrus that are stronger on the left (Ramayya et al. 2010). Cortical folding patterns in the parietal operculum might also be indirectly influenced by asymmetries in the arcuate fasciculus, which links temporal and prefrontal cortex and courses under the parietal operculum (Glasser and Rilling 2008).

We found a leftward asymmetry in PT of only 13% in surface area, which is smaller than some previous reports (Fleschig 1908; Falzi et al. 1982; Kulynych et al. 1993; Habib et al. 1995; Steinmetz 1996). Our cohort included some left handers, which were not analyzed separately, and their inclusion may have reduced the asymmetry. On the other hand, some reports may have overestimated PT surface area asymmetry owing to measurement bias arising from the asymmetry in the 3D orientation of the Sylvian Fissure. A much broader swath of small but significant asymmetries includes left-dominant cortex in insula, anterior temporal and medial parietal and posterior cingulate cortex, and right-dominant cortex in lateral temporal, occipital, and medial prefrontal cortex (Fig. 9E,F). The reproducibility of this pattern needs examination, as does exploration of its possible functional significance. For example, it is of interest to know whether this pattern is correlated with cortical networks revealed in fMRI studies, especially for functions that might be processed differently in the 2 hemispheres.

Anatomical studies of structural asymmetries have reported left dominance of Broca's area gray matter (Falzi et al. 1982) and of probabilistic area 44 (Amunts et al. 1999) as well as a variety of other gray matter asymmetries using voxel-based morphometry (Good et al. 2001; Watkins et al. 2001). However, it is difficult to compare volumetric results with our surface areas asymmetries, given the large methodological differences. The surface-based approach gives a direct quantitative measure that is not affected by statistical differences in the admixtures of gray and white matter at any given volumetric locus.

Cortical Parcellations

Many more human parcellation schemes have been published than were registered to the fs_LR atlas in the present study. Nonetheless, our cross-study comparisons provide several insights. One region of particular interest is area MT, which has been charted using multiple methods. Motion-responsive MT+ has been reported to correspond to hOc5 (Wilms et al. 2005; Amunts et al. 2007) and to an architectonic subdivision delineated using in vivo MRI (Walters et al. 2007). We found substantial overlap between probabilistic retinotopic MT (n = 4) registered via Landmark-SBR and probabilistic hOc5 (n = 10) registered via Energy-SBR. However, determining the precise relationship between these parcels is complicated by differential biases associated with the 2 registration methods. Also, retinotopic MT is smaller than hOc5, but this in part reflects incomplete retinotopic mapping. Malikovic et al. (2007) reported dorsal and ventral architectonic subdivisions (hOc5d and h0c5v) in some hemispheres but not consistently enough to warrant treating them as distinct areas. Altogether, it seems plausible that MT corresponds to one hOc5 subdivision (most likely hOc5d) and that the other architectonic subdivision includes a separate visual area (most likely retinotopic V4t or MSTv).

In general, the combination of spatial uncertainty in identifying areal boundaries, biological variability in the size and location of cortical areas, and variability in registration to an atlas makes it very difficult to equate parcels delineated in different studies unless both parcels are represented by probabilistic maps that have been registered to an atlas using similar if not identical SBR methods. When such maps are available, their boundaries can be estimated objectively by computing the spatial gradient of each probabilistic map and determining whether gradient peaks for each parcel are aligned. This type of approach has been successfully applied in comparing in vivo myelin maps to cytoarchitectonic maps (Glasser and Van Essen 2011) or to gradients in functional connectivity (Glasser et al. 2011).

Our estimate of 150–200 distinct areas per hemisphere in humans exceeds the estimate of ∼140 areas in the macaque reported in the companion study (Van Essen et al. 2011) but is consistent with other estimates for humans (Kaas, 2008). Total surface area per hemisphere is 9.2-fold greater in human (973 cm², present study) versus the macaque (105 cm², Van Essen et al. 2011). The regions of rapid evolutionary expansion in the human lineage lie in lateral temporal, parietal, and prefrontal cortex (Orban et al. 2004; Van Essen and Dierker 2007) and overlap extensively with the unparcellated regions in the human composite map (Fig. 14A). Did rapid evolutionary expansion occur mainly by a differential increase in the size of individual areas in the human lineage, with only a modest increase in the total number of areas? Alternatively, did it occur by the emergence of many new cortical areas, with little change or possibly even a decrease in their average size? Resolving this issue would shed important insights on human brain evolution. However, the problem is especially challenging owing to the inherent difficulty in distinguishing whether heterogeneity in cortical features in a given region reflects 1) small but distinct areas, 2) distinct modules within a larger cortical area, or 3) more irregular fluctuations and gradients within a larger area (Van Essen et al. 2011). In regions of rapid evolutionary change, it is possible that the underlying organizational principles may differ from that in other regions. For example, some regions (e.g., lateral temporal, parietal, an/or prefrontal cortex) might include relatively large cortical areas that contain internal modules or gradients whose arrangement is even more variable than that of evolutionarily more stable regions.

Despite these challenges, better methods of data acquisition and analysis are likely to accelerate progress in addressing issues of cortical parcellation and evolutionary divergence. In vivo neuroimaging methods for systematically and quantitatively charting, structural and functional connectivity offer great potential not only for improving intersubject registration (as noted above) but also for improved cortical parcellation. One general strategy is to identify candidate areal boundaries using sharp transitions in structural and/or functional connectivity based on spatial gradients (Cohen et al. 2008; Nelson et al. 2010) or cluster analysis (Mars et al. 2011). Confidence in such parcellations increases when supported by data from multiple modalities (Nelson et al. 2010; Glasser and Van Essen 2011; Glasser et al. 2011; Mars et al. 2011).

Funding

This work was supported by National Institute of Mental Health (grant R01 MH 60974) and in part by the Human Connectome Project (1U54MH091657-01) from the 16 National Institutes of Health (NIH) Institutes and Centers that Support the NIH Blueprint for Neuroscience Research and by the McDonnell Center for Systems Neuroscience at Washington University.

We thank S. Danker for assistance in manuscript preparation, E. Reid for assistance with data analysis, and Drs J. L. Price, H. Kolster, G. Orban, A. Brewer, B. Wandell, D. Tsao, J. Swisher, S. Pitzalis, and H. Burton for providing data sets used to map cortical parcellations to atlas surfaces. Conflict of Interest: None declared.

References