The human cerebral cortex is made up of a mosaic of structural areas, frequently referred to as Brodmann areas (BAs). Despite the widespread use of cortical folding patterns to perform ad hoc estimations of the locations of the BAs, little is understood regarding 1) how variable the position of a given BA is with respect to the folds, 2) whether the location of some BAs is more variable than others, and 3) whether the variability is related to the level of a BA in a putative cortical hierarchy. We use whole-brain histology of 10 postmortem human brains and surface-based analysis to test how well the folds predict the locations of the BAs. We show that higher order cortical areas exhibit more variability than primary and secondary areas and that the folds are much better predictors of the BAs than had been previously thought. These results further highlight the significance of cortical folding patterns and suggest a common mechanism for the development of the folds and the cytoarchitectonic fields.
The human cerebral cortex is a ribbon of gray matter that is highly folded in order to enable a large surface area to fit in the limited volume provided by the human skull. The folds are intriguing in both their variability and regularity, but little is understood about their relationship to the microstructural organization of the cortex. The cortex itself can be parcellated into a mosaic of microscopically, (i.e., architectonically definable areas) based on localizable and more or less pronounced changes in the laminar distribution of neuronal cell bodies (cytoarchitecture) and/or intracortical myelinated (myeloarchitecture) fibers (Brodmann 1909; Vogt 1911; von Economo 1929; Sarkissov et al. 1955). The most famous of these parcellations is the one proposed by Korbinian Brodmann (Brodmann 1909) a century ago. Most current imaging studies of the human cortex report the location of effects as a “Brodmann area” (BA). This determination is typically made by visual comparison of the functional imaging results with Brodmann's schematic drawings and thus comes with no defined estimate of precision or uncertainty.
Cyto- and myeloarchitectonic differences between adjacent areas that are the basis of the definition of borders between the BAs vary considerably in terms of their subtlety. For example, probably the most salient architectural feature of the cortex is the stria of Gennari, a highly myelinated stripe in layer IV present only in the primary visual cortex (BA 17). The stria of Gennari is one of the few architectural features of the cortex that is detectable in vivo using magnetic resonance imaging (MRI) (Clark et al. 1992; Barbier et al. 2002; Walters et al. 2003). Another prominent cytoarchitectonic feature of the cortex is the layer II islands (Ramon y Cajal 1909; No 1933) in entorhinal cortex (EC, BA 28) that give rise to the perforant pathway through which most of the input from neocortical areas travels to the hippocampus. Using ex vivo MRI at ultrahigh field, we have recently succeeded in robustly visualizing these cell-dense regions throughout the extent of EC (Augustinack et al. 2005). Despite these examples, the vast majority of the architectural characteristics that define borders between adjacent cortical areas are not visible at the resolutions that can be achieved by current neuroimaging technologies. Microscopic analysis of histologically stained brain sections, therefore, still remains the most powerful and reliable tool for cortical parcellation and identification of BAs.
Despite the widespread use of cortical folding patterns to perform ad hoc estimations of the locations of the BAs in individuals, little is understood regarding the relationship of the folds to the BAs or whether there is a hierarchy in the predictability of the BAs. The architectonics are of course important as the mosaic of functionally defined regions that are arrayed across the cortical sheet (e.g., Allman and Kaas 1971; Tootell et al. 1983; Felleman and Van Essen 1991; Sereno and Allman 1991) are strongly linked to the underlying anatomy. Here we use whole-brain histology combined with statistically testable parcellation methods for the identification of cortical areas (Amunts et al. 1999; Zilles et al. 2002) and surface-based analysis (Dale et al. 1999; Fischl, Sereno, Dale 1999; Fischl, Sereno, Tootell et al. 1999) to explicitly test how well the folds predict the locations of the areas. We show that the accuracy with which an area can be predicted from folding patterns appears to be related to its level in the putative cortical hierarchy, with primary and secondary sensory areas being well predicted by surrounding folding patterns, and higher level cognitive areas such as Broca's area (BAs 44 and 45) the most variable with respect to the folds. We anticipate that this type of mapping will allow a more accurate assessment of the uncertainty associated with localization of functional or structural properties of the human brain.
Materials and Methods
Histological Processing Methods
Ten human postmortem brains were processed and analyzed using the techniques described in Schormann and Zilles (1998), Amunts et al. (1999), and Zilles et al. (2002). The silver-stained histological sections of an individual brain were aligned to the postmortem MR volume of the same brain using nonlinear warping (Schormann and Zilles 1998) to build an undistorted 3-dimensional histological volume. The basic steps, which have been employed in numerous studies, for example, Zilles et al. 1995; Geyer et al. 1997; Schormann and Zilles 1998; Amunts et al. 2000; Geyer et al. 2000; Rademacher et al. 2002; Amunts et al. 2005, are as follows.
Histological, cell body–stained sections with cortical regions of interest are imaged under a microscope using a motorized scanning stage and a camera. For subsequent cytoarchitectonic analysis, the gray level index (GLI, [Schleicher and Zilles 1990]) is measured as an index of the volume fraction of cell bodies and GLI images are obtained. Dark pixels correspond to a low volume fraction of cell bodies, light pixels to a high one.
The cortex is covered by intensity line profiles that traverse the cortical ribbon from gray/white boundary to pial surface. The shape of each profile reflects the cytoarchitecture (Schleicher and Zilles 1990).
A distance function is computed to determine the degree of similarity of adjacent blocks of line profiles. A high degree of dissimilarity (or low similarity) indicates a substantial change in the laminar profiles and hence in the underlying cytoarchitecture.
Significant maxima in dissimilarity are those for which the location of the maximum does not depend on the block size but remains stable over large block-size intervals.
Surface-Based Analysis Methods
The reconstructed histological volumes were used to generate surface models of the gray/white interface. This was accomplished in several steps. First, a set of “control points” were manually added to the body of the white matter to guide an intensity normalization step that resulted in the white matter across most of the volume being close to a prespecified value (Dale et al. 1999). This volume was then thresholded and manually edited to separate white matter from other tissue classes. The resulting binary segmentation was used to generate topologically correct and geometrically accurate surface models of the cerebral cortex (Dale et al. 1999; Fischl, Sereno, Dale 1999; Fischl et al. 2001) using a freely available suite of tools (http://surfer.nmr.mgh.harvard.edu/fswiki). An example of the results of this procedure together with the locations of the manually selected control points is given in Figure 1, which shows coronal (top), sagittal (middle), and axial (bottom) slices of a typical volume with the reconstructed gray/white surface shown in yellow. Note that small errors in surface positioning, which would be critical, for example, in a study of cortical thickness, are mostly irrelevant in this study in which we are more concerned with the large-scale geometry of the surface models. The 8 labeled BA maps (areas 2, 4a, 4p, 6, 44, 45, 17, and 18) were sampled onto surface models for each hemisphere, and errors in this sampling were manually corrected (e.g., when a label was erroneously assigned to both banks of a sulcus). A morphological close was then performed on each label. A close of a binary label is a dilation, in which each point that is 0 and neighbors a point that is 1 is set to 1, followed by an erosion, in which each point that is 1 and neighbors a 0 is set to 0. The close was used to remove small holes that arise due to sampling artifacts without distorting the boundary of each label.
The 10 left and 10 right hemispheres were morphed into register using a high-dimensional nonlinear morphing technique that aligns cortical folding patterns (Fischl, Sereno, Tootell et al. 1999). Briefly, this technique maps each individual surface model into a spherical space and then represents the geometry of the surfaces as functions on the unit sphere. The registration of the surfaces is accomplished by maximizing the similarity of these spherical functions, while also constraining the mappings to be invertible and to induce only modest amounts of metric distortion. For these datasets, specifically we used 3 sets of geometric features to drive the registration. The first was the mean curvature of the “inflated” surface (the surfaces displayed in Fig. 2). This was necessary to account for the large-scale geometric distortions present in the data. Next we aligned the “average convexity,” which has been shown to be representative of the primary folding patterns (Fischl, Sereno, Dale 1999). Finally, the mean curvature of the gray/white boundary surface was used as the input feature in order to align secondary and tertiary folds where possible. Each of these features in turn was matched to the corresponding feature in our standard in vivo atlas comprised of 40 subjects distributed in age and pathology (10 with mild Alzheimer's disease). Note that no specific optimization was performed for aligning the BAs presented in this report. Rather, a set of parameters that had been determined to be optimal for aligning primary visual cortex (V1) in a separate ex vivo dataset (Hinds OP, Rajendran N, Polimeni JR, Augustinack JC, Wiggins G, Wald LL, Rosas HD, Potthast A, Schwartz EL, Fischl B, unpublished data) were used with no modification.
In order to quantify the accuracy of the alignment of the underlying BAs, the spherical registration was used to transform each of the 8 BAs for each individual into each of the other individual coordinate systems, and a modified Hausdorff distance was computed. (Note that areas 4a, 4p, and 6 were obtained for only 8 of the 10 total subjects. Each of the other areas was present for every subject.) Specifically, for each point on the boundary of each subject's area in the individual subject space, we computed the minimum distance to the boundary of each other subject projected into the individual subject's original white matter surface model and then computed the average of these. The results of this analysis are displayed in Figure 4. The advantage of this procedure is that it provides a measure in millimeter of the uncertainty of localization and is invariant to the size of an area, a well-known problem for other similarity measures such as the Dice or Jaccard coefficient, which compute the degree of overlap of binary labels, a measure that is affected by the size of the label, with larger labels typically evidencing greater overlap than smaller one.
We constructed spatial probability maps for 8 BAs across 10 human brains (both left and right hemispheres) as shown in Figure 2, which displays the average convexity of the in vivo atlas that is used as a common space. These include the primary and secondary visual areas BA 17 and BA 18, respectively; BA 44 and BA 45 (subdivisions of Broca's area); the somatosensory area BA 2; the primary motor areas 4a and 4p; and finally the premotor area BA 6 (note that these last 3 areas were only analyzed in 8 of the 10 datasets). Frequency estimates of the probability that each point was part of each BA were constructed in a surface-based coordinate system by counting the number of times that a label occurred at a given point and dividing by the total number of subjects for each label. Each point in the surface-based coordinate system can then be probed to determine the probability that it is part of any of the set of labeled BAs.
To assess the accuracy of the surface-based results relative to more standard volumetric procedures, we used the publicly available volumetric probability maps (http://www.fz-juelich.de/inb/inb-3//spm_anatomy_toolbox) constructed using a high-dimensional nonlinear fluid warp (Schormann and Zilles 1998). The accuracy of the 2 techniques was quantified by constructing cumulative histograms of the probability for each nonzero voxel (or vertex) in each probability map for each of the 8 areas, as shown in Figure 3. Each bar represents the probability that a point will be at least that accurate. Because the minimum accuracy would be if the label of only one subject occurred at each location, the smallest value on the x-axis is 0.1 (1 subject out of 10). The histograms always achieve their maximum value of 1 at an accuracy of 0.1 indicating that the entire surface/volume is at least this accurate. Subsequent bars then represent the percent of the surface/volume that is at least this accurate. Thus, the bar at 0.2 represents the portion of the data with an accuracy ≥ 0.2. Ideally, if the normalization perfectly aligns the underlying architectonics, these maps will be binary, with ones in the interior of the region and zeros elsewhere, resulting a flat histogram with the rightmost bin (P ≥ 1) containing as many points as the leftmost (P ≥ 0). The level of the histogram in the high-accuracy bins (more overlap across subjects, toward the right in the histograms) then measures the accuracy with which the underlying coordinate system aligns the borders of the BAs. The accuracy of the surface-based alignment in also aligning the architectonics is summarized in the bottom row, which shows the average of the histograms across the 8 areas (left) and the ratio of the surface and volume histograms on the right. For example, the surface-based coordinate system has greater than 7 times more locations of perfect accuracy than the volumetric one and outperforms the volume at every accuracy level. We believe that this type of result does not reflect the details of the volumetric procedure but rather that surface-based techniques use intrinsically more predictive features—cortical folding patterns—which are not available in the volume. Note that the more commonly used linear alignment procedures (12 parameter affine, not shown) have significantly lower accuracy than the fluid warps.
In order to explicitly quantify how well the folding patterns that were used to construct the surface-based coordinate system predict the locations of the various BAs, we computed the average distance between the boundaries of each individual instance of each BA in its native space to every other individual instance of that BA mapped into that subject's coordinate system, as described in the Materials and Methods section. The results of this analysis are shown in Figure 4. This measure allows both an estimate of the absolute accuracy of localization of each BA as well as a means for comparing how well predicted the boundaries of each BA are relative to the others. Note that errors in the surface reconstructions due to the reduced contrast to noise in the underlying images relative to what can be routinely obtained in vivo only strengthens these findings, as this type of error will only artificially increase our estimates of the variability. Examining Figure 4, it is clear that 1) primary and secondary sensory areas are extremely well predicted by the surrounding geometry and 2) there appears to be progression of accuracy, with the level of predictability diminishing as one moves away from areas devoted to processing sensory inputs and into cortical regions implicated in more cognitive domains.
The most widely used coordinate system in neuroimaging is the one developed by Talairach and Tournoux (Talairach et al. 1967; Talairach and Tournoux 1988), which provides stereotaxic maps for inferring the architectonic localization of cortical effects (e.g., functional or structural differences between populations or conditions). Unfortunately, although popular tools exist for estimating BA from Talairach coordinate (Lancaster et al. 1997; Lancaster et al. 2000), this coordinate system has been shown to be a poor predictor of the locations of both primary sensory (Rademacher et al. 1992; Rademacher et al. 1993; Amunts et al. 2000; Geyer et al. 2000; Morosan et al. 2001; Rademacher et al. 2001) as well as higher order cortical areas (Amunts et al. 1999; Amunts et al. 2005). An alternative and even more widespread approach is to make an ad hoc estimation of the BA containing a given cortical effect by visually comparing individual folding patterns with those in Brodmann's drawings. This procedure, however, is also problematic, because Brodmann's maps are schematized drawings, and thus do not reflect a real individual brain with its folding pattern. Further, Brodmann's drawings give no means of assessing the variability of the relationship between the folds and the cytoarchitectonic boundaries.
The variability of the architectonics has been characterized in several studies, particularly the landmark work of Rajkowska and Goldman-Rakic, in which 7 human left hemispheres were analyzed to characterize the variability in areas 9 and 46 (Rajkowska and Goldman-Rakic 1995a; Rajkowska and Goldman-Rakic 1995b), with reconstructions of the lateral portion of the hemispheres carried out in 5 cases. In this study, considerable variability was found in the morphology of frontal sulcal patterns. Further, by overlaying their architectonic maps on the Talairach atlas, Rajkowska and Goldman-Rakic were able to point out the ambiguity in other published results that reported findings in a particular BA (e.g., an effect reported in area 9 could have been in 45 or 46 also). It has not been clear whether the well-documented inaccuracy of the use of the Talairach coordinate system for localizing BAs reflects the true variability of the underlying architectonic areas or if higher dimensional nonlinear coordinate systems based on other types of macroscopically observable features could be used in order to increase the accuracy of the localization of the underlying cyto- and myeloarchitecture.
In this work, we have shown that computational techniques that explicitly drive folding patterns into register across subjects are also surprisingly accurate at aligning histologically defined BAs, despite having no access to the microscopic properties used to define them. This is particularly true in the primary cortical areas we have investigated, with primary visual cortex (BA 17) being the most predictable, exhibiting in the order of 2.7 mm of median variability in the location of its boundary in both hemispheres across all subjects. In fact, the predictability of all the primary motor and sensory areas that we studied, including BA 17, 4a, 4p (anterior and posterior divisions of BA 4 [Geyer and Ledberg 1996]), and 2 (although recent evidence casts some doubt over whether area 2 should be considered primary or not [Zilles et al. 2004; Toga et al. 2006]), was found to be surprisingly good with a mean uncertainty of approximately 3.7 mm in the surface-based coordinate system. This figure was obtained by computing the median uncertainty of each individual area across each subject and then taking the mean of these. In the few “higher order” areas that we analyzed, the variability increased to 7 mm in the left hemisphere for areas 44 and 45 and 9 mm in the right hemisphere, with significant parts of each area overlapping in all subjects. These core areas of 100% overlap indicate that it should be possible to restrict analysis to regions in which a researcher is confident that an effect is indeed within a given BA, although it is important to note that the geometry of the area will play a role in this type of analysis as well. For example, BAs 44 and 45 exhibit more variability than say BA 4a, but the elongated nature of BA 4a would make it difficult to find many functional MRI voxels solely contained within the predicted location of this cortical area.
Several explanations are possible for this apparent hierarchy in the variability of the location of cortical areas. Variability in position may simply relate to the variability of regional folding patterns as, for example, prefrontal regions are more variable geometrically than perirolandic regions or the region around the calcarine. This, however, begs the question of why primary areas occur near primary folds. If cortical folding patterns are reflective of the tension of subcortical and corticocortical axonal projections (Van Essen 1997), then it may be that the variability in the location of a cortical area relates to the degree of heterogeneity in its pattern of connectivity. Thus, primary areas that are connected to relatively few other cortical areas would be less variable than higher order (multimodal “association”) areas, which project to and receive projections from many more disparate brain regions (Pandya et al. 1988). V1, for example, has connectivity mainly limited to the lateral geniculate nucleus of the thalamus and secondary visual cortex (V2) (for review see Sincich and Horton ). Conversely, area 44 receives major projections from secondary somatosensory area S2 and inferior parietal lobule as well as projections from prefrontal and premotor areas (9, 46v, 47/12, 13, 6), cingulate motor cortex, superior temporal sulcus, and rostral insula (Geschwind 1965; Jones and Powell 1970; Pandya and Yeterian 1996). Area 45 receives its main inputs from superior temporal gyrus (higher auditory cortex) and multimodal areas in the superior temporal sulcus, in addition to other prefrontal areas, somatosensory areas 1 and 2, caudal insula, and visual areas of the inferior temporal cortex (Geschwind 1965; Jones and Powell 1970; Pandya and Yeterian 1996). Variability in cortical localization could thus largely reflect the complexity of the underlying patterns of connectivity, as opposed to being directly related to relative location in a hierarchical arrangement of cortical areas.
It is worth noting that the cytoarchitectonic changes that define the borders between adjacent association cortices (such as 44/45) are considerably more subtle than in primary areas, which typically show reasonably sharp transitions in cellular properties between one area and its neighbors (Van Hoesen 1993), making the precise and repeatable localization of higher areas considerably more difficult. In the present observations, cytoarchitectonic maps were used based on a reliable, observer-independent, and statistically testable microscopical technique (Schleicher et al. 1999), which excludes a systematic increase of variability between primary and higher areas due to such methodical reasons. Phylogenetic factors could play a role in the variability of localization, as it has been posited that primary sensory cortices are the most recent to evolve (Sanides 1970), and therefore, evolutionary age could be reflected in degree of variability. This argument is also supported by the fact that the variability in the volume of neocortical areas 44 and 45 greatly exceeds that of the hippocampus (part of archicortex) (Amunts et al. 1999; Amunts et al. 2005). One important cautionary point is that homologies between macaque and human for areas 44 and 45 have not been definitively established (Deacon 2004). It is also possible that ontogenetic factors influence cortical localization. For example, the order of development could play a role with earlier developing areas being less variable than later ones due to a simple propagation of errors. It is known that primary areas myelinate earlier than higher ones (Flechsig 1901), and there is some evidence that they form earlier as well (Flechsig 1920; Brody et al. 1987), although the early myelination of middle temporal area/V5 would then imply that it would have a stable location with respect to surrounding folding patterns, which does not appear to be the case.
Although the variability across areas is intriguing, one striking feature of our results is the stability of the localization of the BAs with respect to the surrounding folding patterns, as might perhaps have been expected given the demonstrated ability of surface-based registration to align structurally and functionally homologous features of the human cortex (Fischl, Sereno, Tootell et al. 1999; Van Essen 2005). This stability may arise from genetic factors, which are likely to play an important role in the location and size of cortical areas. One prominent hypothesis regarding the formation of cortical areas is that the specification of the architectonic regions is present in a protomap in the proliferative ventricular zone in the form of radial columns that guide the formation and migration of cortical neurons during neurodevelopment (Rakic 1988). There is evidence that the protomap exists without the need for sensory input (e.g., Armentano et al. 2007; Cholfin and Rubenstein 2007), although the size and location of the architectonic areas can be modulated by the modification of afferent input (Goldman-Rakic 1980), perhaps contributing to the observed variability in the localization with respect to the surrounding folding patterns. Thus, the protomap may initially specify the location of the cortical areas with respect to one another, with corticocortical and thalamic connectivity then influencing the creation of cortical convolutions and the final position and size of each architectonic field.
An important implication of the current work is that if the size of a cortical area relates to competence of the functional domain in which the area is implicated, then it may be possible to predict performance levels directly from gross morphology. For example, in recent work, Duncan and Boynton (2003) have shown that visual acuity is predicted by the size of the functionally defined primary visual cortex. Given the accuracy with which the borders of V1 appear to be localized by folding patterns alone, visual acuity should be inferable directly from brain structure. Finally, understanding how the underlying cellular characteristics are arranged with respect to the macroscopically visible folding patterns is an important step in understanding how the folds develop and whether they play a computational role in the processing strategies employed by the human brain.
National Center for Research Resources (P41-RR14075, R01 RR16594-01A1, NCRR BIRN Morphometric Project BIRN002, U24 RR021382); National Institute for Biomedical Imaging and Bioengineering (R01 EB001550, R01EB006758); National Institute for Neurological Disorders and Stroke (R01 NS052585-01); Mental Illness and Neuroscience Discovery Institute, which is part of the National Alliance for Medical Image Computing; National Institutes of Health through the NIH Roadmap for Medical Research (grant U54 EB005149). Helmholtz-Association of Research Centres; the DFG; National Institute of Biomedical Imaging and Bioengineering; National Institute of Neurological Disorders and Stroke; National Institute of Mental Health (to K.Z., K.A.).
Information on the National Centers for Biomedical Computing can be obtained from http://nihroadmap.nih.gov/bioinformatics. Thanks to Jon Polimeni, Verne Caviness, Barry Kosofsky, and Ellen Grant for useful suggestions; Helen Barbas for suggesting the “diversity of input” theory of variability; Kilian Pohl for suggesting the error measure; and Randy Buckner for suggesting the application to the Duncan and Boynton data. Volumetric probabilistic maps of the cytoarchitectonic areas can be downloaded at http://www.fz-juelich.de/inb/inb-3/Home. FreeSurfer and surface-based probabilistic maps of the cytoarchitectonics can be downloaded from http://surfer.nmr.mgh.harvard.edu/fswiki. Conflict of Interest: None declared.