Abstract

We present a technique for automatically assigning a neuroanatomical label to each location on a cortical surface model based on probabilistic information estimated from a manually labeled training set. This procedure incorporates both geometric information derived from the cortical model, and neuroanatomical convention, as found in the training set. The result is a complete labeling of cortical sulci and gyri. Examples are given from two different training sets generated using different neuroanatomical conventions, illustrating the flexibility of the algorithm. The technique is shown to be comparable in accuracy to manual labeling.

Introduction

Techniques for labeling geometric features of the cerebral cortex are useful for analyzing a variety of functional and structural neuroimaging data (Rademacher et al., 1992; Caviness et al., 1996; Paus et al., 1996; Seidman et al., 1997; Goldstein et al., 1999, 2001, 2002; Lohmann et al., 1999). In many cases, what is desired is not a labeling of a discrete set of cortical features, but rather the identification of every point in the entire cortex. This type of labeling is known as a parcellation. Unfortunately, despite their potential utility, cortical parcellations are not commonly used in the neuroimaging community due to the difficult and time-consuming nature of the task of manually parcellating the entire cortex from high-resolution MRI images. While many techniques exist for labeling parts of the cortex (Mangin et al., 1995; Thompson et al., 1996b; Valiant et al., 1996; Sandor and Leahy, 1997; Lohmann, 1998; Lohmann and von Cramon, 1998, 2000; Goualher et al., 1999; Rettmann et al., 2002), none of them provides an automated procedure for generating a complete and detailed labeling of the entire cortex that can incorporate prior information as well as cortical geometry into the labeling.

The inclusion of prior information is a critical feature of a cortical parcellation algorithm. The reason why this is the case is that the divisions that are useful in a cortical parcellation scheme can be based on properties of the brain other than cortical geometry. Many disparate pieces of information, such as knowledge of structure–function relationship, and cytoarchitectonic or receptor labeling properties of regions, may be used by a neuroanatomist in generating a cortical parcellation — information that is not directly available to an automated parcellation procedure from magnetic resonance imaging (MRI) data. A trivial example of this is the fact that many unbroken sulci change names as lobar boundaries are crossed. In addition, functional heterogeneity can make it desirable for adjacent areas to be assigned different labels despite the absence of macroscopic cortical features to distinguish them. Conversely, there are situations in which the secondary folding structure of the cortex (i.e. folds within folds) carries information, such as the location of the hand area in primary motor cortex, which is well predicted by a posterior-pointing secondary fold in the precentral gyrus (Boling et al., 1999).

Various approaches have been taken to the problem of labeling of cortical features. For example, Sandor and Leahy (1997) use a manually labeled atlas brain, which is then warped into correspondence with an individual subject’s surface model. The individual subject surface then inherits the labels from the atlas. Similar surface-fitting approaches have been developed by Thompson et al. (1996b), Valiant et al. (1996), Lohmann (1998) and Lohmann and von Cramon (1998, 2000). A graph-based approach is taken in Goualher et al. (1999), in which sulci are represented by vertices in the graph and arcs connecting them represent their relationship to one another. A manual training set is then used in order to generate a semi-automated procedure for assigning neuroanatomical labels to the detected sulci. Another graph-based approach was taken in Mangin et al. (1995), in which a 3-D skeletonization of cortical sulci that is robust to topological errors is generated and used to construct a graph for sulcal labeling. A somewhat different technique was developed by Lohmann (1998). In this approach, sulci are not represented as abstract quantities such as vertices in a graph, but instead are derived using a watershed transform on an MR volume. This work was extended by Rettmann et al. (2002) to the use of surface models. The segmented regions identified by the watershed transform are then manually labeled by a trained neuroanatomist.

Here we present a technique for using manually labeled data as the basis for an automated parcellation procedure, using an intrinsically cortical coordinate system to store prior statistics and class-conditional densities (Fischl et al., 1999a,b). The use of the manual parcellation as a training set allows the technique to incorporate neuroanatomical convention into the parcellation in regions in which geometry alone is not predictive of a parcellation label. The procedure models the parcellation labels as a first order anisotropic non-stationary Markov random field (MRF), allowing it to capture the spatial relationships between parcellation units that are present in the training set. The anisotropy separates label probabilities in the first and second principal curvature directions, in order to encode the increased probability of differing labels in the direction of maximal curvature. This type of modeling allows the procedure to automatically encode information such as ‘precentral gyrus frequently neighbors central sulcus in the direction of high curvature (moving across the sulcus), but not in the direction of low curvature (along the banks of the sulcus)’. This type of parcellation has commonly been generated by a having a trained anatomist or technician manually label some or all of the structures in the cortex, a procedure that can take up to a week for high-resolution images. Here, we use the results of the manual labeling using the validated techniques of the Center for Morphometric Analysis (CMA) (Rademacher et al., 1992; Caviness et al., 1996) to automatically extract the information required for automating the parcellation procedure. In addition, in order to demonstrate that the technique is sufficiently flexible to encode a different set of neuroanatomical conventions, we show comparable results using an entirely different surface-based (SB) parcellation scheme (Destrieux et al., 1998).

Materials and Methods

Local spatial relationships between labeled structures have been encoded by modeling the labeled image using Markov random fields (MRFs) in a variety of image processing contexts, dating back to their introduction (Geman and Geman, 1984). In the MRF approach, the probability of a label at a given position is computed not just in terms of the data and prior probabilities at that position, but also as a function of the labels in a neighborhood around the point in question. In the context of segmenting MR images, isotropic (all directions are equal) and stationary (the probabilities are the same for all spatial locations) MRFs have been used to provide smoothness constraints on a given segmentation (Held et al., 1997; Kapur et al., 1998; Zhang et al., 2001). In this way, the prior probability of a label is computed by examining how likely the label is given the labels of its neighbors, regardless of the direction of the neighbor, or the position within the brain. While this type of approach can obviate the need for prefiltering of the images, it does not provide for the use of information regarding the spatial relationships of classes to one another. Anisotropic MRFs have been used in the context of segmenting MRI volumes in order to disambiguate the problem of labeling subcortical structures (Fischl et al., 2002), work which is extended here to the labeling of the cortical surface. [Note that the derivations in the following sections closely follow those in Fischl et al. (2002).]

Thus, we base the cortical parcellation procedure on a number of pieces of information incorporated into a space-varying classification procedure. That is, class statistics (e.g. means and covariance matrices) are tabulated regionally throughout an atlas space, using a registration procedure that optimally aligns cortical folding patterns (Fischl et al., 1999b). Prior probabilities are computed via a frequency histogram in the atlas space, allowing the calculation of the probability that each parcellation label occurs at every atlas location. Finally, the prior probability of a given spatial arrangement of parcellation labels is incorporated into the final parcellation procedure. These priors are also computed from a training set, examples of which are given in Figure 1, for each point in the atlas by modeling the parcellation as an anisotropic non-stationary Markov random field, resulting in a procedure that is comparable in terms of accuracy to manual labeling.

Problem Statement

The problem of automatically labeling cortical gyri and sulci can be naturally phrased within the framework of Bayesian parameter estimation theory. In this approach, one can relate the probability of a parcellation P given the observed surface model S to the probability p(S|P) of the geometry of the surface model occurring given a certain parcellation, together with the prior probability of the segmentation:

graphic

The primary advantage of the Bayesian approach is that it allows for the explicit incorporation of prior information via the p(P) term in equation (1). In order to render the problem more tractable in the face of the large degree of variability in cortical folding patterns, both the priors on P and the conditional probability of observing the surface given the classification p(S|P) can be expressed within an atlas space, allowing them to vary as a function of position on the cortical surface (hence making them non-stationary). The advantage of using an atlas space is that coordinates in the atlas have more anatomical meaning than the native coordinate system of the image (Bajcsy et al., 1983; Bookstein, 1989; Miller et al., 1993; Gee et al., 1994; Vannier et al., 1994; Christensen et al., 1996; Collins et al., 1996; Ashburner et al., 1997; Collins and Evans, 1997; Woods et al., 1998; Fischl et al., 1999b; Thompson et al., 2000). Classifiers can then be distributed throughout the atlas, allowing each one to focus on only the small number of classes that may occur within the region for which the classifier is responsible. The number of classes that occur within a region of space is then directly related to the accuracy of the atlas coordinate system. That is, P(P(r) = c) will be 0 for all but a few values of c at each atlas location r. In practice, if the classifiers are reasonably dense in the atlas space then the number of classes at each location is typically relatively small. In this way, the intractable problem of classifying each surface node into one of 80 or so labels is decomposed into a set of tractable problems of classifying the nodes in each region of the surface model into only a small number of labels.

The definition of the atlas requires the calculation of a function f(r), which takes native image coordinate as input, and returns the coordinate of the corresponding point in the atlas. For f to be useful in this context, the coordinates it returns should be related to the anatomical location of r. This type of mapping therefore provides the ability to meaningfully relate coordinates across subjects. In the most general case, we wish to maximize the joint probability of both the parcellation P and the atlas function f:

graphic

The terms p(S|P,f) and p(P|f) in equation (2) provide a natural means for incorporating atlas information into the parcellation procedure. The first term encodes the relationship between the class label at each atlas location and the predicted surface geometry. Using the atlas space, we can allow the class statistics to vary as a function of location, allowing the within-class variations in surface geometry, such as secondary and tertiary cortical folds, to be incorporated into the classification procedure. The second term allows the expression of prior information regarding the spatial structure of the parcellation labels. Finally, the term p(f) provides a means for constraining the space of allowable atlas functions (e.g. continuity, differentiability and invertibility). In the following, we will use the atlas function f described in Fischl et al. (1999b), which establishes a spherical surface-based coordinate system for each cortical hemisphere.

Atlas Construction

In general, two different approaches have been taken to the construction of anatomical atlases from neuroimaging data. The first is to use an individual as a template, and either manually or automatically estimate a transformation that aligns each new dataset with the individual template (e.g. Talairach et al., 1967; Talairach and Tournoux, 1988; Van Essen and Drury, 1997; Van Essen et al., 1998). The alternative is to compile a probabilistic atlas based on the anatomy of a large number of subjects (e.g. Collins et al., 1994; Fox et al., 1994; Mazziotta et al., 1995; Thompson et al., 1997, 2000). Each of these approaches has strengths and weaknesses. The former allows one to represent anatomical structures at as fine a scale as the neuroimaging technology allows, but the atlas is then biased by the idiosyncrasies of the individual anatomy chosen as the template. The latter technique resolves this problem by averaging across anatomies, thus only retaining that which is common in the majority of subjects. Nevertheless, the cross-subject averaging removes potentially useful information in the atlas. Here, we preserve the advantages of each technique by using a group of subjects to construct an atlas that retains information about each parcellation unit at each point in space. Given the atlas function f, and a group of N manually parcellated subjects, we first estimate the prior probability of parcellation label c occurring at each atlas location, independent of all other locations:

graphic

The geometry of the observed surface S at each point is summarized by two quantities encoded as a vector G(r) at each point in the surface (note that here r refers to spherical atlas coordinates, and f(r) refers to individual subject spherical coordinates). The first of these is the average convexity as defined in Fischl et al. (1999a), while the second is the mean curvature of the surface (do Carmo, 1976). Additional quantities such as MRI derived variables (e.g. T1 and T2), Gaussian curvature and cortical thickness can easily be incorporated into this framework. The likelihood of observing the surface geometry G given the parcellation label P(r) is modeled as a Gaussian, the parameters of which are computed from the manually labeled training set:

graphic

where Gi are geometric parameters extracted from each of the set of M surfaces for which label c occurs at location f(r) in the corresponding manually labeled surface models Si [i.e. Si(f(r)) = c]. The covariance matrix for class c at location r in the atlas is then given by:

graphic

Information about surface geometry is thus maintained separately for each parcellation label at every location in the atlas, obviating the need to average geometric information across classes. Finally, we also estimate the pairwise probability that parcellation label c2 is the neighbor at ri when parcellation label c1 is the label at r, for ri ϵ N(r), where N(r) is a neighborhood function of r.

graphic

As before, this information is stored separately for each atlas location. It is important to note that the probabilities are stored separately for each pair of classes as well as for each neighborhood location ri. While it may seem that this would lead to combinatorial explosion and intractable memory requirements, in practice the space is sparse as relatively few configurations of parcellation labels actually occur. Two examples of the manual labeling that is the basis for the atlas are given in Figure 1. The spherical atlas space is made up of a pair of super-sampled icosahedra. The priors are stored at essentially the same resolution as the individual surface models, using an icosahedral model with 163  842 vertices at ∼1 mm spacing, while the class conditional densities are stored 2562 vertices at ∼4 mm spacing, reflecting the slower variation of surface geometry for each label.

Parcellation

As mentioned above, the Bayesian approach allows us to incorporate prior information that is necessary for the parcellation procedure. This prior information takes two forms. The first makes use of the global spatial information provided by f and the atlas in order to express the probability that a given parcellation label occurs at a particular location in the atlas, independent of the local surface geometry. The second encodes the local spatial relationship between parcellation labels, allowing information such ‘precentral gyrus can be anterior, superior or inferior to central sulcus, but never posterior to it’ to be automatically detected and incorporated into the segmentation.

Formally, we compute the maximum a posteriori (MAP) estimate of the parcellation P given an input surface geometry G, and the nonlinear spherical transform f. The MAP estimate can be expressed as maximizing p(P|G,f), the probability distribution of the parcellation given the observed surface geometry and the atlas function. Using Bayes rule we can relate this to the product of the probability of observing the surface geometry with the prior probability of a given spatial configuration of labels p(P):

graphic

Assuming the noise at each vertex is independent from all other vertices in the surface model, we can rewrite p(G|P,f) as the product of the distribution at each vertex over the surface:

graphic

Note that in the more general case in which the spatial correlation structure of the noise is constant across space, the equality in equation (8) should be replaced by proportionality (Worsley et al., 1992; Thompson et al., 1996a), which does not effect any of the subsequent derivations. The distribution of the geometric properties of each class at each location in the atlas is modeled as a Gaussian, the mean vector µc(r) and covariance matrix Σc(r) of which are computed using equations 4 and 5. The probability of observing the geometry at G(f(r)) is then expressed as:

graphic

All that remains is to find an expression for the prior probability of a given parcellation P. Here we assume that the spatial distribution of labels can be well approximated by an anisotropic non-stationary Markov random field. This allows one to encode prior information about the relationship between labels as a function of location within the cortex (i.e. non-stationary), as well as with local direction (i.e. anisotropic). Formally, the Markov assumption can be expressed as:

graphic

That is, the prior probability of a label at a given vertex r is only influenced by the labels within some neighborhood of r. The locality restriction imposed by the Markov model permits the probability of the entire parcellation to be written in terms of neighborhood or clique potentials Vc(P) via the Hammersley–Clifford theorem (Besag, 1974). That is, the probability p(P) can be equivalently characterized by a Gibbs distribution:

graphic

where Z is a normalizing constant and will be dropped in the following, and U(P) is an energy function that can be written in the form:

graphic

The clique potentials Vc(P) encode the energy associated with a certain configuration of labels within the cth clique. Choosing Vc(P) to be log(p(P(r)|P(r1), P(r2) … P(rK)), where r is the central vertex of the cth clique, allows us write the probability of the entire parcellation as the product of the probability of the label at each vertex, given its neighborhood:

graphic

Using Bayes rule, we can rewrite this as:

graphic

Equation (14) allows the probability of a given label to be modulated by any configuration of neighboring labels. While this would be extremely useful, it is unfortunately not computationally tractable to implement, as one would need to compute separate prior probabilities for every combination of neighboring labels that occur. Instead, we make the simplifying assumption that only the first order conditional dependence is important. That is, that the dependence of a label on its neighbors can be expressed as the product of the probability given each of the neighbors:

graphic

where again we have explicitly included the dependence on the neighbor location ri to emphasize that the probability densities are maintained separately for each neighbor position in N(r). Using this assumption, we arrive at an expression for the prior probability of the full parcellation:

graphic

Equation (16) allows two types of prior information to be incorporated into the segmentation procedure. The approximate location a surface feature may occupy within the cortex is given by p(P(r)), which is computed and stored in the atlas using equation (3). The local relationship between parcellation labels is encoded in p(P(ri)|P(r),ri) using equation (6). We currently let the neighborhood function N(r) include the two principal directions on the surface at each location in the atlas space. This allows the parcellation procedure to store information separately for the direction of highest curvature from that of the direction of lowest curvature, as the majority of borders occur in the former.

Directly computing the global MAP estimate of P in equation (7) using the Markov model of equation (16) is computationally intractable. Instead, we employ the iterated conditional modes (ICM) algorithm proposed by Besag (Besag, 1986). In this approach, the parcellation is initialized with the MAP estimate assuming p(P(ri)|P(r),ri) is uniform, as no label has yet been assigned to each vertex. The parcellation is then sequentially updated at each location by computing the label P(r), that maximizes the conditional posterior probability p(P(r)|P(ri),G,ri):

graphic

Equation (17) is then iteratively applied until no vertices are changed. Snapshots of the evolution of this procedure are given in Figure 2. As can be seen convergence is rapid, resulting in a procedure that requires ∼3 min on a 1.5 GHz Pentium III. A small amount of and post-processing is carried out on the results of the labeling to relabel small isolated patches of label to the most likely of the neighboring labels.

Results and Discussion

Here we present the results of applying the surface-based parcellation technique to two different parcellation schemes. The first is the volumetric one in use at the MGH CMA as described in (Rademacher et al., 1992; Caviness et al., 1996), which defines ∼58 separate labels (see Appendix for a complete list), while the second is an intrinsically surface-based (SB) parcellation (Destrieux et al., 1998) based on the conventions established in Duvernoy (1991) with 85 separate parcellation units (see Appendix). For the CMA parcellation we first sample the volumetric labeling onto the reconstructed cortical surface of each subject (Dale et al., 1999; Fischl et al., 1999a, 2001; Fischl and Dale, 2000), the subsequent procedures for the two parcellations are identical.

The data used for the CMA parcellation are part of an ongoing study of schizophrenia (Seidman et al., 1997; Goldstein et al., 1999; Seidman et al., 1999; Goldstein et al., 2002). As part of this study, 36 MRI volumes (2 MP-RAGE scans per subject, motion corrected and averaged) have been manually parcellated by trained technicians. (Note that as the study is not yet complete, we are blinded to diagnosis.) Cortical models were reconstructed for each of the subjects using previously presented techniques, including non-rigid surface-based alignment to a previously constructed spherical atlas (Dale and Sereno, 1993; Dale et al., 1999; Fischl et al., 1999a,b; Fischl and Dale, 2000; Fischl et al., 2001). The manual parcellations were then sampled onto the cortical models at the midpoint of the cortical ribbon using the thickness estimates described in Fischl and Dale (2000). For the SB parcellation, cortical surface models were reconstructed for 12 subjects. Surface-based drawing tools that are part of the FreeSurfer software package (http://surfer.nmr.mgh.harvard.edu) were then used to generate a complete parcellation for all 24 cortical hemispheres. In order to validate the accuracy of the automated parcellation algorithm, we employed a jackknife/leave-one-out technique in which one dataset is removed from each training set, an atlas is constructed, and used to parcellate the left-out set. The error is then quantified point-by-point as the percentage of points at which the manual and automated parcellation are in disagreement. Figure 3 presents the results of this study on an average inflated surface (dark gray are sulcal regions, and light gray are gyral ones). Note that the vast majority of the surface is below threshold (i.e. accuracy >75%) for both the CMA parcellation on the left as well as the SB parcellation on the right. With the median accuracy at ∼81% for the left hemisphere and 83% for the right for the CMA parcellation, and 80% (left) and 79% (right) for the SB (the slightly lower accuracy for the SB parcellation is no doubt due to the greater number of parcellation units in it, some of which are quite small). In addition, the low-accuracy regions are clustered at the boundaries between labels, the precise positions of which are somewhat arbitrary in the manual training set.

The distribution of accuracy collapsed across each cortical hemisphere is given in Figure 4 for the CMA parcellation and Figure 5 for the SB parcellation. The top rows display histograms of the percent correct for the left and right hemispheres, while the bottom row is the cumulative histogram of the same data. As can be seen, almost 40% of the surface was correct 100% of the time for both hemispheres for both of the parcellation schemes.

A scatter plot of the accuracy of the individual subjects is given in Figure 6 CMA parcellation (left) and the SB parcellation (right). Interestingly, a single outlier can be seen in the left hemisphere of both parcellation sets. In both cases this is a subject with a rare folding pattern — a split central sulcus, and even in these case, the accuracy is close to 70%. As our database of manually labeled examples grows, we expect these types of rare neuroanatomies to be handled with better accuracy (the training set that was used for this subjects actually contained no examples of this type of folding pattern, as it was of course left out during construction of the classifiers for labeling it).

In order to directly compare the automated and manual parcellations we computed the mean and standard errors of the surface area of a set of cortical features across all the subjects in each dataset, using the two techniques for each of the parcellation schemes. This comparison is given in Figure 7 for the CMA parcellation, and Figure 8 for the SB parcellation. As can be seen, the surface area for a variety of cortical regions as assessed by the manual (light bars) and automated (dark bars) parcellations is essentially indistinguishable from a statistical standpoint, although it is worth pointing out the uniformly smaller error bars of the automated technique, indicating it is a potentially more powerful technique for detecting variations in cortical surface area.

Finally, a visual comparison of the automated and manual parcellation for the subject with the median accuracy in each the dataset is given in Figure 9. Note that many of the small protuberances in the manual CMA parcellation (on the left) are almost certainly errors in the manual labeling. We are currently developing an intrinsically surface-based adaptation of this parcellation that will minimize this type of error. In addition, the automated parcellation can be used to correct errors in the manual labelings, although of course these datasets could not be used for validation of the automated technique after correction.

Conclusion

Labeling of cortical features is a technique that has been shown to have many neuroscience applications dating back to classic studies on functional lateralization (Geschwind and Levitsky, 1968). More recently, Schlaug and colleagues have found differences in musicians and non-musicians in perirolandic cortex (Schlaug et al., 1995), while Boling et al. have shown folding patterns to be predictive of the location of the motor representation of the hand (Boling et al., 1999). Witelson and Kigar demonstrated asymmetries in the Sylvian fissure that relate to handedness (Witelson and Kigar, 1992), while Goldstein and colleagues carried out an extensive study of sexual dimorphism both cortically and subcortically (Goldstein et al., 2001), as well as the disruption of normal sexual dimorphism in schizophrenia (Goldstein et al., 2002). Other applications include assessments of normal anatomical asymmetry in the sylvian fissure (Rubens et al., 1976), planum temporale changes in neurofibromatosis (Billingsley et al., 2002) and autism (Rojas et al., 2002). Schizophrenia researchers have found effects in many cortical regions including the planum temporale (Barta et al., 1995; Hirayasu et al., 2000), the angular gyrus (Niznikiewicz et al., 2000), the hippocampus (Seidman et al., 1997), as well as the middle frontal gyrus, orbitofrontal cortex, the anterior and paracingulate gyri and the supramarginal gyrus (Goldstein et al., 1999).

While manual methods exist for assessing this type of change, the process of manually parcellating an entire high-resolution structural MR volume requires on the order of a week for a trained neuroanatomist or technician. This makes the routine analysis of large patient and control populations untenable. Finally, manual labeling procedures do not generalize well to the use of multidimensional inputs such as cortical thickness measures, intrinsic MRI tissue parameters and cortical geometry, which can be easily incorporated into the automated parcellation procedure.

The automated method described in this paper for assigning a neuroanatomical label to every point in the cortex has been shown to be comparable in terms of accuracy to a previously validated method of manual parcellation. The accurate labeling of a large number of cortical structures is enabled through the use of both global and local spatial information. The global information is encoded by distributing classifiers throughout an atlas in a surface-based coordinate system and maintaining class statistics on a per-class, per-location basis, allowing the classifiers to be robust to variations cortical folding patterns. Local information is incorporated into the classification procedure by modeling the parcellation as a non-stationary anisotropic Markov random field. The introduction of anisotropy and non-stationarity into the parcellation model allows the spatial relationships of parcellation units to one another to be incorporated into the procedure in a principled fashion.

The automatic parcellation procedure can also be used to automatically define regions of interest (ROIs) for use in functional imaging studies. Specifically, this will allow one to generate average time courses by cortical region, or even parts of regions, facilitating for example the comparison of the response of the motor cortex to that of somatosensory cortex, or anterior superior temporal sulcus (STS) to the posterior STS.

While the results presented here are comparable to those achievable by current manual labeling techniques, we intend to extend them in a number of ways. Primarily, we will investigate the relationship between the manual and automated labelings of each brain. In particular, the CMA parcellation technique is defined in the volume using slice-based viewers, and is currently being updated to use surface-based drawing and visualization tools. When complete, we expect the surface-based extension of the CMA parcellation to be significantly more accurate than the labeling based on the volumetric parcellation. As part of this process, we will perform structure-by-structure analysis of variability both across manual labelings, and between manual and automated labelings for each of the manual training sets. This analysis should allow us to identify regions of low reliability, which will require either updating the definition of the parcellation labels in the case that neither the automated nor the manual labelings is reliable, or the inclusion of additional information into the automated parcellation procedure. This may include the relationship of the cortical regions to subcortical structures, as well as additional cortical information (e.g. regional thickness measures, intrinsic tissue parameters).

The automated nature of the methods described here, in contrast with existing manual or semi-automated techniques, allows for their routine application in large-scale studies. Having access to this type of detailed morphometric information for large populations including various disorders as well as a spectrum of normal controls should facilitate the characterization of the anatomical signatures associated with specific disorders. Ultimately, this may provide a more accurate and sensitive tool for early diagnosis of brain disorders.

Appendix

Labels for CMA cortical parcellation

AG angular gyrus

BFsbcmp basal forebrain

CALC intracalcarine cortex

Cga cingulate cortex, anterior part

CGp cingulate cortex, posterior part

CN cuneus

CO central operculum

F1 superior frontal gyrus

F2 middle frontal gyrus

F3o inferior frontal gyrus, pars opercularis

F3t inferior frontal gyrus, pars triangularis

FMC frontomedial cortex

FO frontal operculum

FOC fronto-orbital cortex

FP frontal pole

H1 Heschl’s gyrus

INS insula

LG lingual gyrus

MW medial wall

OF occipital fusiform gyrus

Oli occipital lateral gyri, inferior part

Ols occipital lateral gyri, superior part

OP occipital pole

PAC paracingulate cortex

PCN precuneus

Pha parahippocampal gyrus, anterior part

PHp parahippocampal gyrus, posterior part

PO parietal operculum

POG postcentral gyrus

PP planum polare

PRG precentral gyrus

PT planum temporale

SC subcallosal cortex

SCLC supracalcarine cortex

Sga supramarginal gyrus, anterior part

SGp supramarginal gyrus, posterior part

SMC supplementary motor cortex or juxtapositional lobule cortex

SPL superior parietal lobule

T1a superior temporal gyrus, anterior part

T1p superior temporal gyrus, posterior part

T2a middle temporal gyrus, anterior part

T2p middle temporal gyrus, posterior part

T3a inferior temporal gyrus, anterior part

T3p inferior temporal gyrus, posterior part

Tfa temporal fusiform gyrus, anterior part

TFp temporal fusiform gyrus, posterior part

TO2 middle temporal gyrus, temporo-occipital part

TO3 inferior temporal gyrus, temporo-occipital part

TOF temporal fusiform gyrus, temporo-occipital part

TP temporal pole

UN unknown

Labels for surface-based parcellation

Corpus_callosum

G_and_S_Insula_ONLY_AVERAGE

G_cingulate-Main_part

G_frontal_inf-Opercular_part

G_frontal_inf-Triangular_part

G_frontal_superior

G_insular_long

G_occipital_inferior

G_occipital_superior

G_occipit-temp_med-Lingual_part

G_orbital

G_parietal_inferior-Angular_part

G_parietal_superior

G_precentral

G_rectus

G_subcentral

G_temp_sup-Lateral_aspect

G_temp_sup-Planum_tempolale

G_temporal_middle

Lat_Fissure-ant_sgt-ramus_horizontal

Lat_Fissure-post_sgt

Pole_occipital

S_calcarine

S_central_insula

S_cingulate-Marginalis_part

S_circular_insula_inferior

S_collateral_transverse_ant

S_frontal_inferior

S_frontal_superior

S_intermedius_primus-Jensen

S_intraparietal–and_Parietal_transverse

S_occipital_inferior

S_occipital_superior_and_transversalis

S_occipito-temporal_medial_and_S_Lingual

S_orbital_medial-Or_olfactory

S_paracentral

S_pericallosal

S_precentral-Inferior-part

S_subcentral_ant

S_suborbital

S_supracingulate

S_temporal_superior

S_transverse_frontopolar

Unknown

This research is funded jointly by the National Center for Research Resources (R01 RR16594-01A1), the National Institute of Neurological Disorders and Stroke, the National Institute of Mental Health and the National Cancer Institute (R01-NS39581, R01-NS34189). Further support was provided by the National Center for Research Resources (P41-RR14075 and R01-RR13609), the National Institute of Mental Health (NIMH MH56956), the Theodore and Vada Stanley Foundation, the Amytotrophic Lateral Sclerosis Association, Faculté de Médecine/CHRU, Tours, France, the Mental Illness and Neuroscience Discovery (MIND) Institute and the National Association for Research in Schizophrenia and Depression. Thanks to Steve Hodge and Jennifer Koch for providing much-needed support and organization. Thanks also to Andrea Boehland, Bobbie Holland and Nicole Cullen for providing manual parcellations for use in this paper.

Figure 1. Lateral views of five left hemisphere examples of the CMA (top two rows) and SB (bottom two rows). Each pair of inflated (above) and pial (below) surface representations is from the subject.

Figure 1. Lateral views of five left hemisphere examples of the CMA (top two rows) and SB (bottom two rows). Each pair of inflated (above) and pial (below) surface representations is from the subject.

Figure 2. Snapshots of the ICM algorithm. From left to right: initial configuration, first, second, and final configuration (after 14 iterations).

Figure 2. Snapshots of the ICM algorithm. From left to right: initial configuration, first, second, and final configuration (after 14 iterations).

Figure 3. Map of per cent incorrect for the CMA (left) and the surface-based (right) parcellations. Top: lateral view, middle: medial view, bottom: ventral view. Note that the majority of the surface is correct >75% of the time in both parcellation schemes.

Figure 3. Map of per cent incorrect for the CMA (left) and the surface-based (right) parcellations. Top: lateral view, middle: medial view, bottom: ventral view. Note that the majority of the surface is correct >75% of the time in both parcellation schemes.

Figure 4. Total per cent accuracy for the automated CMA parcellation (integrated across each cortical hemisphere). Top row: histogram of per cent accuracy for the left (left) and right (right) hemispheres. Bottom row: cumulative histograms for the same.

Figure 4. Total per cent accuracy for the automated CMA parcellation (integrated across each cortical hemisphere). Top row: histogram of per cent accuracy for the left (left) and right (right) hemispheres. Bottom row: cumulative histograms for the same.

Figure 5. Total per cent accuracy for the SB parcellation (integrated across each cortical hemisphere). Top row: histogram of percent accuracy for the left (left) and right (right) hemispheres. Bottom row: cumulative histograms for the same.

Figure 5. Total per cent accuracy for the SB parcellation (integrated across each cortical hemisphere). Top row: histogram of percent accuracy for the left (left) and right (right) hemispheres. Bottom row: cumulative histograms for the same.

Figure 6. Scatter plot of accuracy for CMA (left) and surface-based (right) parcellations. The median accuracies are 78 and 79%, respectively.

Figure 6. Scatter plot of accuracy for CMA (left) and surface-based (right) parcellations. The median accuracies are 78 and 79%, respectively.

Figure 7. Comparison of mean and standard error of the surface area of 10 of the CMA parcellation units across 36 subjects (left hemisphere on the left, and right hemisphere on the right) for the manual (white bars) and automated (black bars) parcellations. FP, frontal pole; INS, insula; F1, superior frontal region; F2, middle frontal region; F3t, inferior anterior middle frontal region; F3o, inferior posterior middle frontal region; PRG, precentral gyrus; TP, temporal pole; T1a, anterior superior temporal sulcus; T1p, posterior superior temporal sulcus.

Figure 7. Comparison of mean and standard error of the surface area of 10 of the CMA parcellation units across 36 subjects (left hemisphere on the left, and right hemisphere on the right) for the manual (white bars) and automated (black bars) parcellations. FP, frontal pole; INS, insula; F1, superior frontal region; F2, middle frontal region; F3t, inferior anterior middle frontal region; F3o, inferior posterior middle frontal region; PRG, precentral gyrus; TP, temporal pole; T1a, anterior superior temporal sulcus; T1p, posterior superior temporal sulcus.

Figure 8. Comparison of mean and standard error of the surface area of 10 of the surface-based parcellation units across all 12 subjects (left hemisphere on the left, and right hemisphere on the right) for the manual (white bars) and automated (black bars) parcellations: Scal, calcarine sulcus; Scen, central sulcus; Spost, postcentral sulcus; Gpre, precentral gyrus; Gcun, cuneus; Sst, superior temporal sulcus; GsF, superior frontal gyrus; Spo, parieto-occipital sulcus; Scing, cingulate sulcus; SsF, superior frontal sulcus.

Figure 8. Comparison of mean and standard error of the surface area of 10 of the surface-based parcellation units across all 12 subjects (left hemisphere on the left, and right hemisphere on the right) for the manual (white bars) and automated (black bars) parcellations: Scal, calcarine sulcus; Scen, central sulcus; Spost, postcentral sulcus; Gpre, precentral gyrus; Gcun, cuneus; Sst, superior temporal sulcus; GsF, superior frontal gyrus; Spo, parieto-occipital sulcus; Scing, cingulate sulcus; SsF, superior frontal sulcus.

Figure 9. Manual (left) and automated (right) parcellation results for the subject with the median accuracy. Top, middle and bottom rows are lateral, medial and ventral views, respectively.

Figure 9. Manual (left) and automated (right) parcellation results for the subject with the median accuracy. Top, middle and bottom rows are lateral, medial and ventral views, respectively.

References

Ashburner J, Neelin P, Collins DL, Evans AC, Friston KJ (
1997
) Incorporating prior knowledge into image registration.
Neuroimage
 
6
:
344
–352.
Bajcsy R, Lieberson R, Reivich M (
1983
) A computerized system for the elastic matching of deformed radiographic images to idealized atlas images.
J Comput Assist Tomogr
 
5
:
618
–625.
Barta P, Pearlson G, McGilchrist I, Lewis ATRW, Pulver A, Vaughn D, Casanova GMF, Powers R (
1995
) Reversal of asymmetry of the planum temporale in schizophrenia.
Am J Psychiatry
 
152
:
715
–721.
Besag J (
1974
) Spatial interaction and the statistical analysis of lattice systems (with discussion).
J R Stat Soc Appl Stat Ser B
 
36
:
192
–326.
Besag J (
1986
) On the statistical analysis of dirty pictures (with discussion).
J R Stat Soc Ser B
 
48
:
259
–302.
Billingsley RL, Schrimsher GW, Jackson EF, Slopis JM, Moore BD III (
2002
) Significance of planum temporale and planum parietale morphologic features in neurofibromatosis type 1.
Arch Neurol
 
59
:
616
–622.
Boling W, Olivier A, Bittar RG, Reutens D (
1999
) Localization of hand motor activation in Broca’s pli de passage moyen.
J Neurosurg
 
91
:
1
–12.
Bookstein FL (
1989
) Principal warps: thin-plate splines and the decomposition of deformations.
IEEE Trans Pattern Anal Machine Intell
 
11
:
567
–585.
Caviness VSJ, Meyer J, Makris N, Kennedy DN (
1996
) MRI-based topographic Parcellation of the human neocortex: an anatomically specified method with estimate of reliability.
J Cogn Neurosci
 
8
:
566
–587.
Christensen GE, Rabbitt RD, Miller MI (
1996
) Deformable templates using large deformation kinematics.
IEEE Trans Image Process
 
5
:
1435
–1447.
Collins D, Evans A (
1997
) Animal: validation and applications of non-linear registration-based segmentation.
IJPRAI
 
11
:
1271
–1294.
Collins DL, Neelin P, Peters TM, Evans AC (
1994
) Data in standardized Talairach space.
J Comput Assist Tomogr
 
18
:
292
–205.
Collins DL, Le Boualher G, Benugopal R, Caramanos A, Evans AC, Barillot C (
1996
) Cortical constraints for non-linear cortical registration. In: Visualization in biomedical computing (Hohne KH, Kikinis R, eds), pp.
307
–316. Berlin: Springer.
Dale AM, Sereno MI (
1993
) Improved localization of cortical activity by combining EEG and MEG with MRI cortical surface reconstruction: a linear approach.
J Cogn Neurosci
 
5
:
162
–176.
Dale AM, Fischl B, Sereno MI (
1999
) Cortical surface-based analysis I: segmentation and surface reconstruction.
Neuroimage
 
9
:
179
–194.
Destrieux CE, Halgren E, Dale AM, Fischl B, Sereno MI (
1998
) Variability of the human brain studied on the flattened cortical surface.
Soc Neurosci Abstr
 
24
:1164.
do Carmo M (
1976
) Differential geometry of curves and surfaces. Englewood Cliffs, NJ: Prentice-Hall.
Duvernoy H (
1991
) The human brain. Vienna: Springer-Verlag.
Fischl B, Dale AM (
2000
) Measuring the thickness of the human cerebral cortex from magnetic resonance images.
Proc Natl Acad Sci USA
 
97
:
11044
–11049.
Fischl B, Liu A, Dale AM (
2001
) Automated manifold surgery: constructing geometrically accurate and topologically correct models of the human cerebral cortex.
IEEE Trans Med Imaging
 
20
:
70
–80.
Fischl B, Salat D, Busa E, Albert M, Dieterich M, Haselgrove C, van der Kouwe A, Killiany R, Kennedy D, Klaveness S, et al. (
2002
) Whole brain segmentation. Automated labeling of neuroanatomical structures in the human brain.
Neuron
 
33
:
341
–355.
Fischl B, Sereno MI, Dale AM (
1999
) Cortical surface-based analysis ii: inflation, flattening, a surface-based coordinate system.
Neuroimage
 
9
:
195
–207.
Fischl B, Sereno MI, Tootell RBH, Dale AM (
1999
) High-resolution inter-subject averaging and a coordinate system for the cortical surface.
Hum Brain Mapp
 
8
:
272
–284.
Fox PT, Mikiten S, Davis G, Lancaster JL (
1994
) BrainMap: a database of human functional brain mapping. In: Functional neuroimaging (Thatcher RW, Hallett M, Zeffiro T, John ER, Huerta M, eds), pp.
95
–106. San Diego: Academic Press.
Gee JC, Haynor DR, Reivich M, Bajcsy R (
1994
) Finite element approach to warping of brain images. In: Proceedings of SPIE Medical Imaging 1994: Image processing (Loew MH, ed.). Bellingham: SPIE.
Geman S, Geman D (
1984
) Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images.
IEEE Trans Pattern Anal Machine Intell
 
6
:
721
–741.
Geschwind N, Levitsky W (
1968
) Human brain: left–right asymmetries in temporal speech region.
Science
 
161
:
186
–187.
Goldstein JM, Goodman JM, Seidman LJ, Kennedy DN, Makris N, Lee H, Tourville J, Caviness VS Jr, Faraone SV, Tsuang MT (
1999
) Cortical abnormalities in schizophrenia identified by structural magnetic resonance imaging.
Arch Gen Psychiatry
 
56
:
537
–547.
Goldstein JM, Seidman LJ, Horton NJ, Makris N, Kennedy DN, Caviness VS Jr, Faraone SV, Tsuang MT (
2001
) Normal sexual dimorphism of the adult human brain assessed by in vivo magnetic resonance imaging.
Cereb Cortex
 
11
:
490
–497.
Goldstein JM, Seidman LJ, O’Brien LM, Horton NJ, Kennedy DN, Makris N, Caviness VS Jr, Faraone SV, Tsuang MT (
2002
) Impact of normal sexual dimorphisms on sex differences in structural brain abnormalities in schizophrenia assessed by magnetic resonance imaging.
Arch Gen Psychiatry
 
59
:
154
–164.
Goualher GZ, Procyk E, Collins DL, Venugopal R, Barillot C, Evans AC (
1999
) Automated extraction and variability analysis of sulcal neuroanatomy.
IEEE Trans Med Imaging
 
18
:
206
–217.
Held K, Kops ER, Krause BJ, Wells WM III, Kikinis R, Muller-Gartner H-W (
1997
) Markov random field segmentation of brain mr images.
IEEE Trans Med Imaging
 
16
:
876
–886.
Hirayasu Y, McCarley RW, Salisbury DF, Tanaka S, Kwon JS, Frumin M, Snyderman D, Yurgelun-Todd D, Kikinis R, Jolesz FA, et al. (
2000
) Planum temporale and Heschl gyrus volume reduction in schizophrenia — a magnetic resonance imaging study of first-episode patients.
Arch Gen Psychiatry
 
57
:
692
–699.
Kapur T, Grimson WEL, Wells WM, Kikinis R (
1998
) Enhanced spatial priors for segmentation of magnetic resonance imagery. MICCAI’98. 1st International Conference on Medical Image Computing and Computer Assisted Intervention, MIT, Cambridge, MA, October 11–13.
Lohmann G (
1998
) Extracting line representations of sulcal gyral patterns in MR images of the human brain.
IEEE Trans Med Imaging
 
17
:
1040
–1048.
Lohmann G, von Cramon DY (
1998
) Sulcal basins and sulcal strings as new concepts for describing the human cortical topography. IEEE Workshop on BioMedical Image Analysis, June 27–27, Santa Barbara, CA, p. 24.
Lohmann G, von Cramon DY (
2000
) Automatic labelling of the human cortical surface using sulcal basins.
Med Image Anal
 
4
:
179
–188.
Lohmann G, von Cramon DY, Steinmetz H (
1999
) Sulcal variability of twins.
Cereb Cortex
 
9
:
754
–763.
Mangin J-F, Frouin V, Bloch I, Regis J, Lopez-Krahe J (
1995
) From 3-D magnetic resonance images to structured representations of the cortex topography using topology preserving deformations.
J Math Imag Vision
 
5
:
297
–318.
Mazziotta JC, Toga AW, Evans AC, Fox P, Lancaster J (
1995
) A probabilistic atlas of the human brain: theory and rationale for its development.
Neuroimage
 
2
:
297
–318.
Miller MI, Christensen GE, Amit Y, Grenander U (
1993
) A mathematical textbook of deformable neuroanatomies.
Proc Natl Acad Sci USA
 
90
:
11944
–11948.
Niznikiewicz M, Donnino R, McCarley RW, Nestor PG, Iosifescu DV, O’Donnell B, Levitt J, Shenton ME (
2000
) Abnormal angular gyrus asymmetry in schizophrenia.
Am J Psychiatry
 
157
:
428
–436.
Paus T, Otaky N, Caramanos Z, MacDonald D, Zijdenbos A, D’avirro D,
Gut
 tmans D, Holmes C, Tomaiuolo F, Evans AC (
1996
) In vivo morphometry of the intrasulcal gray matter in the human cingulate, paracingulate, and superior-rostral sulci: hemispheric asymmetries, gender differences and probability maps.
J Comp Neurol
 
376
:
664
–673.
Rademacher J, Galaburda AM, Kennedy DN, Filipek PA, Caviness VS (
1992
) Human cerebral cortex: localization, parcellation and morphometry with magnetic resonance imaging.
J Cogn Neurosci
 
4
:
352
–374.
Rettmann ME, Han X, Xu C, Prince JL (
2002
) Automated sulcal segmentation using watersheds on the cortical surface.
Neuroimage
 
15
:
329
–344.
Rojas DC, Bawn S, Benkers TL, Reite ML, Rogers SJ (
2002
) Smaller left hemisphere planum temporale in adults with autistic disorder.
Neurosci Lett
 
328
:
237
–240.
Rubens AB, Mahowald MW, Hutton JT (
1976
) Asymmetry of lateral (sylvian) fissures in man.
Neurology
 
26
:
620
–624.
Sandor S, Leahy R (
1997
) Surface-based labeling of cortical anatomy using a deformable atlas.
IEEE Trans Med Imaging
 
16
:
41
–54.
Schlaug G, Jäncke L, Huang YX, Steinmetz H (
1995
) In vivo evidence of structural brain asymmetry in musicians.
Science
 
267
:
699
–701.
Seidman LJ, Faraone SV, Goldstein JM, Goodman JM, Kremen WS, Matsuda G, Hoge EA, Kennedy D, Makris N, Caviness VS, et al. (
1997
) Reduced subcortical brain volumes in nonpsychotic siblings of schizophrenic patients: a pilot magnetic resonance imaging study.
Am J Med Genet
 
74
:
507
–514.
Seidman LJ, Faraone SV, Goldstein JM, Goodman JM, Kremen WS, Toomey R, Tourville J, Kennedy D, Makris N, Caviness VS, et al. (
1999
) Thalamic and amygdala-hippocampal volume reductions in first-degree relatives of patients with schizophrenia: an MRI-based morphometric analysis.
Biol Psychiatry
 
46
:
941
–954.
Talairach J, Tournoux P (
1988
) Co-planar stereotaxic atlas of the human brain. New York: Thieme Medical Publishers.
Talairach J, Szikla G, Tournoux P, Prosalentis A, Bordas-Ferrier M, Covello L, Iacob M, Mempel E (
1967
) Atlas d’anatomie stereotaxique du telencephale. Paris: Masson.
Thompson P, Schwartz C, Lin RT, Khan AA, Toga AW (
1996
) Three-dimensional statistical analysis of sulcal variability in the human brain.
J Neurosci
 
16
:
4261
–4274.
Thompson PM, Schwartz C, Toga AW (
1996
) High-resolution random mesh algorithms for creating a probabilistic 3-D surface atlas of the human brain.
Neuroimage
 
3
:
19
–34.
Thompson PM, MacDonald D, Mega MS, Holmes CJ, Evans AC, Toga AW (
1997
) Detection and mapping of abnormal brain structure with a probabilistic atlas of cortical surfaces.
J Comput Assist Tomogr
 
21
:
567
–581.
Thompson P, Woods R, Mega M, Toga A (
2000
) Mathematical/computational challenges in creating deformable and probabilistic atlases of the human brain.
Hum Brain Mapp
 
9
:
81
–92.
Valiant M, Davatzikos C, Bryan BN (
1996
) Finding 3-D parametric representations of the deep cortical folds. IEEE Mathematical Methods in Biomedical Image Analysis workshop, San Francisco, CA, June 21–22.
Van Essen DC, Drury HA (
1997
) Structural and functional analyses of human cerebral cortex using a surface-based atlas.
J Neurosci
 
17
:
7079
–7102.
Van Essen DC, Drury HA, Joshi S, Miller MI (
1998
) Functional and structural mapping of human cerebral cortex: solutions are in the surfaces.
Proc Natl Acad Sci USA
 
95
:
788
–795.
Vannier MW, Miller MI, Grenander U (
1994
) Modeling and data structures for registration of a brain atlas of mutimodality images. In: Functional neuroimaging — technical foundations (Thatcher RW, Hallett M, John ER, Huerta E, eds), pp.
217
–221. New York: Academic Press.
Witelson SF, Kigar DL (
1992
) Sylvian fissure morphology and asymmetry in men and women: bilateral differences in relation to handedness in men.
J Comp Neurol
 
323
:
326
–340.
Woods RP, Grafton ST, Holmes CJ, Cherry SR, Mazziotta JC (
1998
) Automated image registration: II., Inter-subject validation of linear and nonlinear models.
J Comput Assist Tomogr
 
22
:
155
–165.
Worsley KJ, Evans AC, Marrett S, Neelin P (
1992
) A three dimensional statistical analysis for CBF activation studies in human brain.
J Cereb Blood Flow Metab
 
12
:
900
–918.
Zhang Y, Brady M, Smith S (
2001
) Segmentation of brain MR images through hidden Markov random field model and the expectation maximization algorithm.
IEEE Trans Med Imaging
 
20
:
45
–57.