Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Hum Brain Mapp. Author manuscript; available in PMC 2010 March 31.
Published in final edited form as:
PMCID: PMC2847686

Region-of-Interest-Based Analysis With Application of Cortical Thickness Variation of Left Planum Temporale in Schizophrenia and Psychotic Bipolar Disorder


In neuroimaging studies, spatial normalization and multivariate testing are central problems in characterizing group variation of functions (e.g., cortical thickness, curvature, functional response) in an atlas coordinate system across clinical populations. We present a region-of-interest (ROI)-based analysis framework for detecting such a group variation. This framework includes two main techniques: ROI-based registration via large deformation diffeomorphic metric surface mapping and a multivariate testing using a Gaussian random field (GRF) model on the cortical surface constructed by the eigenfunctions of the Laplace–Beltrami operator. We compared our GRF statistical model with a pointwise hypothesis testing approach, whose P-value is corrected using false discovery rate or random field theory at several smoothness scales. As an illustration, we applied this framework to a clinical study of the cortical thickness of the left planum temporale (PT) in subjects with psychotic bipolar disorder, schizophrenia, and healthy comparison controls. Our results show that the anterior portion of the left PT is thinner in the psychotic bipolar and schizophrenic groups than in the healthy control group, and the posterior portion of the left PT shows the reversal finding. Moreover, there may be a greater thickness variation in the left PT in psychotic bipolar patients when compared with that in schizophrenic patients.

Keywords: LDDMM, Laplace–Beltrami operator, Gaussian random field, cortical thickness, planum temporale, schizophrenia, psychotic bipolar disorder


Unlike other diseases (e.g., heart disease), mental illnesses cannot be definitively identified by physiology—for example, through a blood test. Diagnosis of mental illnesses is often made on the basis of symptoms, course of illness, and family history when available. Several mental diseases or disorders have psychotic symptoms, which do not distinguish these diseases. For instance, patients with bipolar disorder often experience hallucinations and delusions so that they may be incorrectly diagnosed as having schizophrenia. One aim of neuroimaging studies is to quantify brain anatomy and function that can serve as complementary information to potentially help diagnose mental diseases. These measurements could be scalar, such as gray matter (GM) volume, cortical surface area, or function indexed over anatomical coordinates (e.g., image volume, cortical surface), such as cortical thickness, functional response, and curvature. Compared to scalar measurement, the function indexed over the anatomical coordinates provide rich information that allows to characterize the variation of the function at points on the anatomical structure across clinical populations, which makes it attractive in current neuroimaging studies. However, it requires studying anatomical variation across subjects and performing statistical analysis on high-dimensional data relative to a small number of populations. All of these require two main areas of technical developments: spatial normalization of the anatomies and multivariate statistical testing.

Most existing spatial normalization methods have focused on registration of the whole brain structure, either for image volumes or hemispherical cortical surfaces. These are powerful methods and provide first insight of disease effects on the brain structure and functions, without requiring prior knowledge of which specific subregions are related to dysfunctions. However, it is a challenge to improve the quality of spatial normalization due to the complicated nature of the brain structure and answer how function alteration is distributed in an affected subregion of the cortex. Several studies have reported that surface-based representation of functions is superior to volume-based representation for studying subregions of the cortex [Van Essen et al., 1998], since the cortical surface preserves the convoluted nature of the brain gyral–sulcal pattern. Moreover, region-of-interest (ROI)-based spatial normalization increases the statistical power compared to the whole brain mapping approach, as shown in several functional magnetic resonance imaging (fMRI) studies [Miller et al., 2005; Nieto-Castanon et al., 2003], and high-dimensional mapping approaches provide more accurate spatial normalization than low-dimensional mapping approaches. This suggests that the high-dimensional spatial normalization methods are needed to register ROI-based cortical surfaces.

Besides spatial normalization, characterization of functions indexed over the curved cortical surface is another challenge in hypothesis testing of group difference across clinical populations. Existing approaches are to first perform hypothesis testing at each point on the cortical surface (e.g. t-test) and then determine the threshold using the random field theory (RFT) [Chung et al., 2005; Worsley, 1994; Worsley et al., 1996, 2004] or false discovery rate (FDR) [Genovese et al., 2002] to correct P-value for multiple comparisons. The statistical results from these approaches are significantly influenced by the smoothness of data.

We introduce a new framework of ROI-based analysis on the functions indexed over a subregion of the cortical surface M [subset or is implied by] R3 in clinical populations. We present the ROI-based surface normalization method via large deformation diffeomorphic metric mapping (LDDMM), which is a high-dimensional spatial normalization approach with properties of one to one, invertible smooth with smooth inverse, and topologically preserving maps. The LDDMM-surface mapping is a fully automatic registration method, where the surface is represented as a geometric object using its normal vectors. In the LDDMM framework [Dupuis et al., 1998; Joshi and Miller, 2000; Trouvé, 1995], the diffeomorphic map from one surface to an atlas is computed. We then characterize observed functions F(·) on the ROI M via a Gaussian random field (GRF) in the atlas. Unlike the approach using the RFT for the correction of multiple comparisons [Chung et al., 2005; Worsley, 1994; Worsley et al., 1996, 2004], we directly construct a GRF model and describe it through its orthonormal representation. Thus, function F(x), x [set membership] M can be written as F(x) = ΣkFkψk(x), where ψk are the eigenfunctions of the Laplace–Beltrami (LB) operator that encodes the geometric structure of the cortex (e.g., area, angle, and length). Fk is the Gaussian random variable associated with ψk. Here, we do not assume that the representation F1, F2, … is exact (i.e., that F(·) can be reconstructed from it), but we expect that it will conserve the necessary information to perform statistical analysis, which will amount to performing inference on F(x) reduced to the coefficients Fi, i=1, 2, …. Therefore, a direct benefit from the GRF orthonormal representation is the reduction of dimensionality of the data. The smoothness of data is also automatically determined by the number of the LB eigenfunctions used in the statistical analysis. We compared our GRF with the pointwise hypothesis testing approach using FDR or RFT correction for multiple comparisons at several smoothness scales.

As an illustration, we apply this analysis framework to detect thickness variation in schizophrenia and psychotic bipolar disorder on the left planum temporale (PT). The left PT is selected because it is the auditory association cortex responsible for language and speech processing function and has received great attention in schizophrenia research [Hirayasu et al., 2000; Kasai et al., 2003; Kwon et al., 1999; Shapleske et al., 2001]. This clinical study aims to investigate the similarity and distinction of thickness distribution over the left PT in schizophrenia and psychotic bipolar disorder relative to healthy comparison controls.


Figure 1 illustrates the flow chart for data processing and statistical analysis. In the following, we describe each step in this flow chart using the left PT structure as an example.

Figure 1
Flow charts of data processing.

Subjects and Image Acquisition

Sixty subjects were selected from schizophrenia and bipolar disorder studies at the Division of Psychiatric Neuroimaging at Johns Hopkins University School of Medicine. The subjects were assigned to three groups for neuroanatomical comparison, subjects with schizophrenia (20), psychotic bipolar disorder (20), as well as healthy comparison subjects (20). All subjects with right-handness [Annett, 1970] were selected into the three groups in terms of age and gender match (Table I). The diagnosis of each individual was determined by the consensus of a research psychiatrist who conducted a semistructured interview and a research assistant who used the criteria from DIGS/MINI [Sheehan et al., 1998] or SCAN/CIDI-SF [Kessler et al., 1998]. Patients with bipolar disorder were particularly selected in the sense that they also have a positive history of psychotic symptoms, such as hallucination, delusion, and thought disorder. Given the size of the sample, we did not consider the effects of medication or examine details of the symptom profiles in our analysis. All subjects gave informed consent for their participation after the risks and benefits of participation were explained to them prior to MRI scanning.

Mean and standard deviation of age within each group

The subjects were examined by means of high-resolution magnetic resonance (MR) 1.5 T Philip scans, which were acquired by using MPRAGE sequence (repetition time = 13.40 ms, echo time = 4.6 ms, flip angle = 20, number of acquisition = 1, matrix 256 × 256), with 1-mm3 isotropic resolution across the entire cranium. Raw MR data were reformatted from signed 16-bit to unsigned 8-bit using ANALYZE [Robb et al., 1989].

Data Processing

Segmentation and surface generation

The skull was removed from each MRI volume using watershed module in ANALYZE [Robb et al., 1989]. A 3D ROI subvolume (see contour region in Fig. 2a) encompassing the superior temporal gyrus (STG) was automatically masked for the left hemisphere in each subject, after manually defining three landmarks in an MRI volume (anterior and posterior commissures, midpoint between them) [Ratnanather et al., 2003]. Bayesian segmentation using the expectation–maximization algorithm to fit the compartmental statistics was used to label voxels in the subvolume as gray matter (GM), white matter (WM), or cerebrospinal fluid (CSF) [Joshi et al., 1999; Miller et al., 2000]. Then surfaces were generated at the GM/WM interface using a topology-correction method and a connectivity-consistent isosurface algorithm [Han et al., 2001, 2002]. The surface location in the MRI volume is illustrated by cyan line in Figure 2a, and the cortical surface is shown in Figure 2b. In Figure 2b, the boundary of PT was delineated by tracking principal curves from the retroinsular end of the Heschl’s Sulcus (HS, red line) to STG (blue line), along the posterior STG up to the start of the posterior ramus and back to the retroinsular end of the HS (green line) via dynamic programming [Honeycutt et al., 2000; Ratnanather et al., 2003]. We cut off the posterior ramus at the final (steepest) upward ascent of the ascending ramus (i.e., where the ascending ramus begins) [Honeycutt et al., 2000; Ratnanather et al., 2003]. Three landmarks on each surface (intersections between STG and HS, between STG and posterior ramus, between posterior ramus and HS) were manually defined in this procedure and formed the triangulated shape of PT (Fig. 2c). Each PT surface was represented by a triangulated mesh with ~1,000 vertices (Fig. 2c). This procedure consistently delineates the PT structure across subjects and has been validated [Ratnanather et al., 2003].

Figure 2
(a) One sagittal slice of a MRI image. The cyan curve is the gray and white matter boundary. (b) Cortical surface containing the superior temporal gyrus (STG). (c) Planum temporale (PT) surface with three boundaries (STG by blue, anterior boundary by ...

PT volume extraction

To extract the PT subvolume, we computed the shortest Euclidean distance of voxel in the ROI volume to the surface M shown in Figure 2b given by

equation M1

where i indexes the ith voxel in the ROI volume and vj is a vertex of surface M. We augmented voxel i with distance di and label li (the index of vertex that voxel i is the closest to). Based on the information of li, the volume associated with the PT surface shown in Figure 2c was extracted from the ROI volume. Its distance map is shown in Figure 2d. The red line is the location of the PT surface. The region of GM and CSF is colored by green while the WM region is encoded by blue.

Cortical thickness maps

We studied thickness via local labeled cortical mantle distance map (LLCMD), which is a modification of LCMD developed by Miller et al. [2003]. The LLCMD maps π(vi) can be studied from the distribution on distances of GM voxels locally around each vertex vi on surface M to approximate the cortical thickness for this local coordinate. Assume Lj to be label of tissue types (WM, GM, or CSF) for each voxel j. Ni is a set of vertices that contains neighbors of vertex vi. Vi = {j|lj = vk, ∀k [set membership] Ni [union or logical sum] {i}}, a set of voxels labeled as closest to vertex vk, where k is a neighbor of vi. The LLCMD map for vertex vi is defined as the frequency of occurrence of GM voxels belonging to vi and its neighbors as a function of cortical location d given by

equation M2

where 1 is an indicator function. Then, the cortical thickness at vertex vi is defined as the 95 percentile of π(vi) in order to remove the error caused by the tissue misclassification at the significance level α = 0.05. In summary, the cortical thickness map is a function of vertex location vi on the cortical surface.


Our variational problem for registering two cortical surfaces is posed in the LDDMM framework in the form of

equation M3

where the first term is a regularization term to guarantee [var phi]t, the flow of template object Itemp at time t, is diffeomorphic (one to one, smooth with smooth inverse map) and its velocity, vt, is an element of Hilbert space of smooth vector fields, V. The template object flows from the identical map (id) and stops at the closest shape to the target object, Itarg. The second term is a matching functional, D([var phi]1 · Itemp, Itarg), that quantifies the closeness between the deformed object [var phi]1 · Itemp and target object Itarg at time t = 1.

To briefly describe LDDMM-surface mapping algorithm under the LDDMM framework, we assume the cortical surface embedded in R3 to be a two-dimensional manifold, in the sense that the neighborhood of every point on the surface is equivalent to a two-dimensional plane in Euclidean space. Such a plane can be uniquely defined by a point and a vector originated at this point and normal to the plane. We define a matching functional, D([var phi]1 · Itemp, Itarg), for registering surfaces in diffeomorphic setting based on their position and normal vectors. Let Itemp = {cf, ηf} and Itarg = {ch, ηh} be the template and target triangulated meshes represented by center points of triangles on the surface and their corresponding normal vectors. Denote the deformed template surface [var phi]1 · Itemp = {cf1, ηf1}, where equation M4 is the center of deformed triangle f and equation M5 is the normal vector to deformed triangle f at location cf1. Let f, g be indices of triangles on the surface Itemp and h, q be indices of triangles on the surface Itarg. The data attachment term D([var phi]1 · Itemp, Itarg) in Eq. (1), which quantifies the closeness between [var phi]1 · Itemp and Itarg, is given in the form

equation M6

where kW (x, y) is a kernel and in practice defined as an isotropic Gaussian kernel matrix, equation M7. ||xy|| is Euclidean distance between points x and y and id3×3 is a 3 × 3 identical matrix. The first two terms are intrinsic energies of the two surfaces [var phi]1 · Itemp and Itarg. The last term gives penalty to mismatching between normal vectors of Itarg and those of [var phi]1 · Itemp [Vaillant and Glaunès, 2005; Vaillant et al., 2007]. From the definition of D([var phi]1 · Itemp, Itarg), this mapping approach does not require any predefined correspondence information (e.g., paired landmarks or curves) and works for both open and closed surfaces.

To investigate diagnostic group difference in thickness, we transformed all PT surfaces into a template using LDDMM-surface mapping. One left PT surface, shown in Figure 4a, from the control group was selected as the template, because its shape has a typical PT triangle structure and its surface area is close to the average surface area in the population. We then carried the thickness maps to the template, based on the information from the closest distance of vertices on the template to vertices on targets.

Figure 4
The surface matching results. Panel (a) shows the template PT surface. Panels (b–g) illustrate the template PT and its deformed PT structures overlayed with the original PT. Red, cyan, and green represent the template, original, and deformed template ...

GRF Model on the Cortical Surface

We constructed a GRF model indexed over the cortical surface for the purpose of performing statistical inference on the cortical thickness. Assume the cortical surface to be an oriented two-dimensional manifold M in R3. Let H(M) be the Hilbert space of square integrable and real-valued function on M with inner product

equation M8

where dv(x) is area measure at the location of x [set membership] M. Using the inner product operation in H(M), we can define the real-valued random field F(x) on M as a GRF with mean μF(x) [set membership] H(M) and covariance operator KF(x, y) if for any u(x) [set membership] H(M), left angle bracket u, F right angle bracket is normally distributed with mean μu = left angle bracket u, μF right angle bracket and variance equation M9. The GRF is completely characterized by its mean μF and covariance operator KF. Assume that the GRF is quadratic mean continuous everywhere ( equation M10) with continuous covariance operator in H(M). Then the orthogonal representation of the GRF in H(M) is given by

equation M11

where ψk is a real-valued orthonormal basis function, spanning the space of square integrable fields on H(M). Fk, k = 1, 2,… are independent, Gaussian random variables with fixed means E[Fk = μk and variances var(Fk) = σ2, where equation M12 and equation M13. Therefore, the mean and covariance operator of F(x) are given as the forms associated with the span {ψk(x)}:

equation M14

equation M15

The covariance operator KF is often empirically defined from observed data. Its principal components from principal component analysis serve as orthonormal basis functions in H(M), which are data dependent [Joshi et al., 1997]. To avoid this problem, we directly derive a complete set of orthonormal basis from a self-adjoint linear differential operator L on H(M), such that its eigenfunctions construct a complete set of orthonormal basis functions in the space of H(M). We particularly chose L as the LB operator on the cortical surface, which is the extension of the Laplacian from the regular grid to a surface. It is interpreted by the intrinsic properties of the surface (e.g., angle, length, and area). We computed its eigenfunctions via the LB spectral problem with Neumann boundary conditions for surface M posed as

equation M16

Δ is the LB operator. u is the local coordinate with entries u1 and u2 on M. n is the normal vector on the boundary of M. Let λ and ψ(u) be the LB eigenvalue and its associated eigenfunction. Define E(ψ(·)) as

equation M17

Based on the divergence theorem and variational principle, it can be shown that the solutions of ψ in (5) are critical points of E(ψ(·)), and λ is its stationary values. Therefore, instead of solving the partial differential equation in (5), we directly solved the variational problem in (6) via the finite element method and no special triangulation procedure was required [Qiu et al., 2006]. In this setting, the eigenvalue problem of the LB operator becomes the algebraic eigenvalue problem. We sorted λ in an increasing order such that λ0 ≤ λ1 ≤ λ2 ≤ …. The numerical solution for ψ0 is trivial due to the fact that λ0 = 0 and gives the null space of the LB operator. The computational details were described elsewhere [Qiu et al., 2006]. As an example, Figure 3 shows ψ1− ψ6 eigenfuctions on a left PT cortical surface. As the order of the eigenfunction increases, the positive and negative valued regions alternate rapidly, indicating the high-frequency signal. The LB eigenfunctions are relatively template–independent, in the sense that the LB eigenfunctions computed from different PTs appear similarly up to rotation due to the multiplicity of eigenvalues [Qiu, 2006].

Figure 3
The first six orthonormal bases of the Laplace–Beltrami operator for the left PT surface in panels (a)–(f).

Gaussian random variables Fk in (2) are determined by the inner product of function F(x) with basis ψk(x) in the form

equation M18

where F(x) is a function indexed over the extrinsic template coordinates.

Statistical Analysis on the Cortical Surface

We developed a statistical analysis method to investigate nonuniform change in thickness on the cortical surface across clinical populations. In the following, F(x) refers to the cortical thickness map on the cortical surface after sub-stracting the mean. Assume H(M) is a Hilbert space of GRF with covariance operator defined as equation M19, where M represents the template surface and ψk(x) is the kth LB eigenfunction on M. From the previous section, a GRF F(x) [set membership] H(M) was described in the form of (2). GRFs {f1i, f2i, i = 1, 2, … m} are expanded as a linear combination of eigenfunctions ψk (x):

equation M20

where the first letter in the superscript represents a group index while the second is a subject index. m is the total number of subjects in each group and n represents the number of coefficients we study. Assume equation M21 and equation M22 be the vectors of the expansion coefficients. Assume that the {f1i, f2i, i = 1, 2, …, m} are samples from GRFs, then the expansion coefficients equation M23 are samples from Gaussian distributions with means { equation M24} and the common sample variance Σ. The hypothesis under the earlier assumptions equation M25 against the alternative equation M26, which detects nonuniform thickness change between two groups distributed on the cortical surface. Hotelling’s T2 is chosen as statistic in this multivariate hypothesis testing, which is a generalization of Student’s t statistic. If Hotelling’s T2 is greater than a threshold at a significance level of 0.05, then the null hypothesis is rejected. We did not perform the correction for multiple comparisons, partly because the degree of freedom of the Hoteling’s T2 test counts the size of Fn. Under the null hypothesis H0 it is shown that

equation M27

which determines P-values [Anderson, 1958].


Examples of Surface Matching

Figure 4 shows several examples before and after the LDDMM-surface mapping. The original PT surfaces in panels (b,d,f) colored by cyan are also shown on the left column of Figure 5, while the original PT surfaces in panels (c,e,g) are illustrated on the right column of Figure 5. Red and green surfaces are, respectively, the template and its deformed versions on each panel. The results suggest that the deformed template is much closer to the original PT shape compared to one before matching. The accuracy of the mapping procedure has been quantitatively evaluated on the same dataset using Euclidean surface to surface distance and geometric measurements as well as the comparison with the landmark matching [Vaillant et al., 2007]. In the Euclidean positional validation, after the LDDMM-surface mapping, about 90% of vertices on the template have distances to all other PTs less than the MRI resolution of 1 mm. We reported that the diffeomorphic maps from the LDDMM-surface mapping algorithm carries more of the variability of the anatomical structures compared with the landmark mapping algorithm.

Figure 5
Examples of cortical thickness maps on left planum temporales. From the top to bottom rows, thickness maps are examples from control (CON), schizophrenia (SCZ), and bipolar (BPD) groups, respectively.

Cortical Thickness Maps

Figure 5 illustrates examples of the cortical thickness maps on left PT surfaces for subjects chosen from healthy control, schizophrenic, and psychotic bipolar groups. Cortical thickness maps of left PT indicate that thickness varies across the PT from 1.5 to 4 mm. The PT is thin at the bottom of HS, then progressively thicker away from HS, and finally thinner toward the posterior ramus. The results suggest that the variation in cortical thickness may be related to the geometry of the PT cortical surface, in the sense that regions with high curvature (perhaps gyral regions) have larger thickness value while regions with low curvature (sulcal regions) have smaller thickness value. Our estimated thickness is consistent with the findings from postmortem studies [Chance et al., 2004; von Economo, 1929]. Thickness is the smallest at the anterior boundary (the bottom of HS, roughly about 1.5 mm) and then it progressively increases in the direction parallel to the anterior boundary toward to STG.

To investigate variation in thickness across subjects, a natural way is to first study the correspondence between anatomical structures using transformation. In this study, we deformed the template PT to all other left PT surfaces and then remapped the cortical thickness maps to the template. The average thickness and standard deviation maps within each group are illustrated in Figure 6. The average thickness maps in control, schizophrenic, and psychotic bipolar groups show the same trend as we observed in individuals. Sulcal regions are thinner and gyral regions are thicker, which is consistent with the cytoarchitechtonic study from von Economo [1929]. The standard deviation maps suggest that the variation of thickness over the PT surface is around 0.25–0.50 mm in these three groups. However, the variation in sulcal region (near HS) is high in the healthy control group, which may partly due to variation of sulcal anatomy.

Figure 6
Panels (a–c), respectively, show average thickness maps within control, schizophrenia, and bipolar groups. Panels (d–f) illustrate standard deviation of thickness maps in control, schizophrenia, and bipolar groups.

Pairwise Comparisons

Figure 7a–c show pairwise differences in the average thickness maps between control and schizophrenia, control and psychotic bipolar, as well as schizophrenia and psychotic bipolar disorder. Red denotes the region thicker in the former group than in the latter group while blue represents the region thinner in the former group than in the latter group. For instance, the red region in panel (a) indicates that the thickness is larger in the healthy comparison control group than in the schizophrenia group. The blue region indicates that the thickness is smaller in the control than in the schizophrenia group.

Figure 7
Panels (a–c), respectively, show pairwise differences of control vs. schizophrenia, control vs. bipolar, and schizophrenia vs. bipolar. For instance, red color denotes the region that thickness is greater in the control group than in the schizophrenic ...

Scalar comparison

In order to detect thickness change over all left PT structures due to effects of schizophrenia and psychotic bipolar disorder, we first reduced the thickness map, a function of location on the cortical surface, into a scalar (mean thickness over the surface). The boxplots in Figure 8 shows median values in three groups (red lines) are close to each other even though it shows the order that the median value of the control group is the greatest and one of the schizophrenic group is the smallest. We performed one-way ANOVA for comparing the means in the three groups of the scalars. The statistical P-value (P = 0.5976) suggests that all three groups exhibit no statistically significant group differences in mean thickness over the surface. This implies that psychotic bipolar disorder and schizophrenia have little effect on the average thickness over the left PT. Assume that the GM ribbon is a laminar structure with surface area S and average thickness T. Then the GM volume, V, can be approximated by the product of S and T, V = S × T. No average thickness reduction with the GM volume reduction from previous neuroimaging morphological studies thus implies that the left PT surface area is reduced by rule of equation M28 in both the schizophrenia and psychotic bipolar groups [Hirayasu et al., 2000; Kasai et al., 2003; Kwon et al., 1999]. The surface area reduction in the patient groups may indicate increased minicolumnar density (the number of minicolumns, elementary functional units normal to the cortical surface, per unit surface area) in schizophrenia [Buxhoevenden et al., 2000], which results in reduced minicolumnar segregation (the surface area per minicolumn) in the PT. This has been proposed to be related to functional asymmetry in language processing. Loss of segregation would limit the effective resolution of the language processing and reduce its discriminative capacity, such as reduced differentiation between semantic categories and between close and distant word associations [Chen et al., 1994; Weisbrod et al., 1998].

Figure 8
Boxplots of mean thickness values within control (CON), schizophrenia (SCZ), and bipolar (BPD). Red line is the location of median value within each group. Bottom and top blue lines are lower and upper quartiles. *Represents individual measurement.

Surface-based comparison

In a further study, we investigated nonuniform effects on thickness due to the occurrence of psychotic bipolar disorder and schizophrenia on the left PT cortical surface. To do so, the Hotelling’s T2 test was performed on coefficient vectors as described previously. The null hypothesis is that the coefficient vectors have no difference between two groups, which implies that the group difference in thickness does not statistically significantly contain the component of their corresponding LB eigenfunctions. Table II lists T2 statistics and the corresponding P-values. We progressively increased the number of coefficients during each test upto five coefficients (approximated by square root of the number of subjects in one group). These P-values were determined from F-distribution. The P-values suggest that the statistical significance in distinct pattern of the thickness map between healthy comparison control and schizophrenia groups is caused by the LB eigenfunction shown in Figure 3a. Such a pattern can be reconstructed by subtracting the mean coefficient in the schizophrenia group from one in the healthy comparison control, and then multiplying this eigenfunction (Fig. 7d). Similarly, the pattern of statistically significant difference in the thickness map between the healthy comparison control and psychotic bipolar groups is due to the three eigenfunctions shown in Figure 3a–c. Its reconstruction is shown in Figure 7e using these three eigenfunctions. Our results show the trend that the posterior regions of left PT is thicker and the anterior region of left PT (near HS) is thinner in both psychotic bipolar and schizophrenic groups, compared to healthy comparison subjects. This suggests that the GM loss in the left PT most occurs in the anterior portion between STG and HS that may be related to sulcal widening during the disease development.

The T2 test statistics and the corresponding P-values for the tests involving the first n coefficients

We do not see statistically significant evidence of thickness difference between schizophrenia and psychotic bipolar groups even though the paired comparisons of the patient groups with the healthy control group would suggest that the difference in thickness between the groups of psychotic bipolar disorder and schizophrenia would exhibit the pattern similar to the LB eigenfunctions shown in Figure 3b,c if there was any. The similar pattern of the thickness alternation between the two patient groups may be partly due to the fact that all patients with bipolar disorder in this study have psychotic symptoms (hallucination, delusion, and thought disorders) as part of their illness, which are also common psychotic symptoms in schizophrenic patients. The subtle altered thickness is in agreement with a relative preservation of normal cytoarchitecture of the left PT in the schizophrenia and psychotic bipolar groups from the previous postmortem study [Beasley et al., 2005]. That is in contrast to findings of reduced neuronal size, glial cell density, and increased neuronal density in prefrontal and anterior cingulate cortices in these disorders [Cotter et al., 2001, 2002, 2004], even though there is evidence for decreased neuronal clustering in the left PT in both schizophrenia and bipolar disorder [Cotter et al., 2004]. This suggests a subtle deviation of neuronal arrangement rather than density and size, which may underlie the thickness alteration pattern and functional abnormalities previously observed in the superior temporal cortex in schizophrenia and bipolar disorder. However, the reconstructed thickness difference maps in Figure 7d,e suggest that there may be a greater thickness variation of the left PT in psychotic bipolar patients when compared with those in the schizophrenic group.

Comparison With Pointwise Based Approach

Pointwise hypothesis testing has been widely used in neuroimaging studies, including volume-based and surface-based fMRI, thickness, GM density, and so on [e.g. Chung et al., 2005; Genovese et al., 2002]. FDR [Genovese et al., 2002] and the RFT [Chung et al., 2005; Worsley et al., 1996, 2004] approaches are used to find objective and effective threshold for pointwise statistics derived from neuroimaging data. However, the threshold determined by both approaches is dependent on the smoothness of data. In this section, we first explore thickness difference between the healthy comparison control and schizophrenia groups from the pointwise hypothesis testing approach thresholded using FDR or RFT at different smooth scales. Then, we compare the findings with those using our method.

In the pointwise hypothesis testing, we first applied the heat kernel smoothing procedure [Chung et al., 2005] to our thickness data on the left PT template surfaces using smoothing kernel of full width half maximum (FWHM) at levels of 20 and 30 mm. Rows (a,b) in Figure 9 illustrate the average thickness maps over the healthy comparison control and schizophrenia groups at different smooth levels. As FWHM increases, the thickness value at each location of the left PT template surface gets close to each other. However, the heat kernel smoothing procedure may not guarantee to converge to the sample mean of the thickness data on a surface with boundary when a large FWHM is applied [Chung et al., 2005]. We next hypothesized that the thickness at each location of the left PT template is equal in both healthy comparison control and schizophrenia groups against the alternative that the thickness is larger in the healthy comparison control group than in the schizophrenia, based on the previous neuroimaging finding that the PT reduces the GM volume in schizophrenia [Hirayasu et al., 2000; Kasai et al., 2003; Kwon et al., 1999]. We performed two-sample t-test on the thickness data at each vertex of the left PT template surface. The t-value threshold was obtained using both FDR at a desired FDR bound of 0.10 and the RFT at a significance level of 0.05 when considering smoothness of the data in multiple comparisons (Table III). For instance, when FWHM is chosen as 20 mm, the t-value threshold is 2.08 via FDR and 2.40 via RFT. Vertices where the t-value is above the threshold represent the locations where the thickness is larger in the healthy comparison control group than in the schizophrenia group. Row (c,d) in Figure 9 respectively show the t-value maps thresholded using FDR or RFT. The t-values below the threshold were replaced by zero and remained otherwise. Red denotes the region where the thickness is larger in the healthy control group than in the schizophrenia group. Each column corresponds to one level of the data smoothness. Both FDR and RFT gave roughly the same region with statistically significantly larger thickness in the control group than in the schizophrenia group. However, as the smoothness increases, such a region becomes smaller via FDR (see row (c)) while it becomes larger via RFT (see row (d)). This observation is consistent with the previous report [Genovese et al., 2002] that suggested the FDR method becomes more conservative as correlations increase while the RFT method is in contrast typically more conservative for unsmoothed data. In summary, the statistical conclusion drawn from the FDR or RFT methods is dependent on the smoothness of the data. As a comparison with our method, the smoothness of the data is automatically taken into account in the GRF model and determined by the number of eigenfunctions in the model. Our statistically significant result shown in Figure 7d captures the pattern that is observed in Figure 7a and shows that thickness is larger in the control than in the schizophrenia in the same region as FDR and RFT detected as well as some other region on the PT template surface. In this sense our method increases the statistical power.

Figure 9
Rows (a,b) show average thickness maps within the healthy comparison control and schizophrenia groups at different smooth scales (FWHM = 20, 30 mm) via heat kernel smoothing [Chung et al., 2005]. Row (c) shows the t-value maps where t-value remains when ...
The t-value thresholds determined by FDR with a desired FDR bound at 0.10 or RFT at a significance level of 0.05 at several levels of the smoothness


This paper presents a novel framework of ROI-based analysis on functions indexed over the surface of cortical subregions, such as cortical thickness, curvature, and functional response. We applied this framework to examine the cortical thickness variation of the left PT in schizophrenia and psychotic bipolar disorder. This framework has three major advantages. First, ROI-based spatial normalization is superior to whole brain mapping approach and increases the statistical power, which has been demonstrated in several functional MRI studies [Miller et al., 2005; Nieto-Castanon et al., 2003]. The LDDMM-surface mapping used in this study provides higher quality of registration than other ROI-based spatial normalization methods (e.g., landmark or curve mappings). More discussion has been reported in its validation paper [Vaillant et al., 2007]. Second, orthonormal data representation via the LB eigenfunctions reduces the data dimensionality to circumvent the need of the correction for multiple comparisons in pointwise hypothesis testing approaches and to examine global effects of the data. Third, our method concludes the finding that cortical thinning is variably distributed across the PT cortical surface in subjects with schizophrenia. It may help to shed light on the variable results that have been reported by other investigators. Notably, some groups have reported GM volume reduction or thinning in the PT [Hirayasu et al., 2000; Kasai et al., 2003; Kwon et al., 1999], while others have failed to find such differences [Shapleske et al., 2001]. To the extent that the cortical thickness of the PT variable in its anatomical location across patients with schizophrenia and psychotic bipolar disorder, small differences in defining the anatomical boundaries of the PT could easily lead to differing results [see details in Honeycutt et al. 2000].

Compared to the spherical brain mapping approaches [Fischl et al., 1999; Thompson et al., 2004; Van Essen, 2004, 2005; Van Essen and Drury, 1997; Van Essen et al., 2001; Tosun et al., 2004a,b], LDDMM-surface does not require an intermediate spherical representation of the brain. This intermediate step would introduce large distortion of the brain structure, which does not appear consistently across subjects. As a consequence, mappings would begin with such a distortion error. However, our approach directly works on cortical surfaces and does not require surfaces with correct topology. Furthermore, our matching approach can map two open surfaces from one to the other and the boundaries of two surfaces need not match if the geometries near the boundary are quite different from one another. This is particularly attractive because it allows us to study variation of anatomies due to effects of diseases in cortical substructures.

The LB eigenfunctions of an arbitrary surface extend the previous work on surface harmonic representation for shapes such as curved manifolds. Earliest work in computer vision was that of Pentland and Sclaroff on surface harmonics in the deterministic setting [Pentland and Sclaroff, 1991]. Joshi et al. [1997] then extended this for stochastic representations such as appropriate for morphometric inference being pursued by the Csernansky group in the hippocampus [Csernansky et al., 1998, 2004]. The empirical basis functions were computed as principal components for the cortical surface and then statistical testing was performed on coefficients associated with these principal components [Joshi et al., 1997]. The natural question becomes whether the principal components would differ after removing or adding data into analysis. If it is the case, then statistical conclusions strongly depend on the estimated components. Furthermore, the principal components may present the variation caused by mismatching anatomical structures when the matching procedure does not provide accurate results, since the components are computed in data-driven fashion. However, the LB eigenfunctions with Neumann boundary conditions overcome these drawbacks. First, it does not rely on data but only on geometry of the cortical surface. Second, the first few eigenfunctions capture slowly varying pattern of data over the cortical surface but not the variation due to mismatching, since the mismatch error (Figs. 4 and 5d in Vaillant et al. [2007]) is randomly distributed over the surface, which corresponds to higher order eigenfunctions. Therefore, the LB eigenfunctions not only can be used as a smoothing basis on the cortical surface, which is equivalent to the Fourier basis on the regular grid [Qiu et al., 2006], but also can serve as a good representor for detecting nonuniform anatomical change. Furthermore, the GRF model constructed via the LB eigenfunctions is superior in the sense that the statistical testing does not require the correction for multiple comparisons and the data smoothness is automatically taken into account by the LB eigenfunctions, which are major concerns in current brain statistical analysis packages, such as statistical parametric mapping [Chung et al., 2005; Worsley, 1994; Worsley et al., 1996, 2004]. Moreover, for a surface with boundary, the kernel smoothing [Chung et al., 2005] does not guarantee the convergence in the sense that average thickness over the surface is constant when applying different smoothing kernel sizes. Thus, it makes difficult to estimate the smoothness level of the data when the cortical surface with boundary is considered as ROI (e.g. PT). However, the limitation of our GRF model is that it requires the computation of the LB eigenfunctions. For a large surface, the computation is intensive. One way to overcome this issue is to apply iterative methods for eigenvalue problems of large sparse symmetric matrices in order to speed up computation [Golub and van der Vorst, 2000]. Furthermore, the LB eigenfunction is a global representation of signals indexed over the cortical surface, and a set of basis functions with local support is more useful in the hypothesis testing. However, computing a set of basis functions with local support only dependent on the surface geometry is difficult and unclear yet and needs to be further investigated.

In this paper, our interest has only focused on investigating the thickness alteration of the left PT in schizophrenia and psychotic bipolar disorder. Compared to the right PT, the left PT is dominant at processing speech, particularly at processing shorter temporal transitions in the speech signal [Harasty et al., 2003; Seldon, 2005; Zatorre et al., 2002]. Furthermore, in word comprehension, the left hemisphere is associated with activation of narrow semantic fields and the right hemisphere with wide semantic fields [Nakagawa, 1991]. The question is whether these functional differing capabilities between the left and right hemispheres are reflected in the thickness of this PT region in healthy subjects and psychotic patients. To go further, we must examine asymmetry of thickness in patients compared to controls. Therefore, future studies on these questions are still needed. Certainly, the presented ROI-based analysis framework provides the possibility to do so.


Contract grant sponsor: NIH; Contract grant numbers: R01 MH064838, R01 EB00975, P41 RR15241.


  • Anderson TW. An Introduction to Multivariate Statistical Analysis. New York: Wiley; 1958.
  • Annett M. A classification of hand preference by association analysis. Brit J Psychol. 1970;61:303–321. [PubMed]
  • Beasley C, Chana G, Honavar M, Landau S, Everall IP, Cotter D. Evidence for altered neuronal organisation within the planum temporale in major psychiatric disorders. Schizophr Res. 2005;73:69–78. [PubMed]
  • Buxhoevenden D, Roy E, Switala A, Casanova MF. Reduced interneuronal space in schizophrenia. Biol Psychiatry. 2000;47:681–683. [PubMed]
  • Chance SA, Tzotzoli PM, Vitelli A, Esiri MM, Crow TJ. The cytoarchitecture of sulcal folding in Heschl’s sulcus and the temporal cortex in the normal brain and schizophrenia: Lamina thickness and cell density. Neurosci Lett. 2004;367:384–388. [PubMed]
  • Chen EYH, Wilkins AJ, Mckenna PJ. Semantic memory is both impaired and anomalous in schizophrenia. Psychol Med. 1994;24:193–202. [PubMed]
  • Chung MK, Robbins SM, Dalton KM, Davidson RJ, Alexander AL, Evans AC. Cortical thickness analysis in autism with heat kernel smoothing. NeuroImage. 2005;25:1256–1265. [PubMed]
  • Cotter D, Mackay D, Landau S, Kerwin R, Everall IP. Reduced glial cell density and neuronal size in the anterior cingulate cortex in major depressive disorder. Arch Gen Psychiatry. 2001;58:545–553. [PubMed]
  • Cotter D, Mackay D, Chana G, Beasley C, Landau S, Everall IP. Reduced neuronal size and glial cell density in area 9 of the dorsolateral prefrontal cortex in subjects with major depressive disorder. Cereb Cortex. 2002;12:386–394. [PubMed]
  • Cotter D, Mackay D, Frangou S, Hudson L, Landau S. Cell density and cortical thickness in Heschl’s gyrus in schizophrenia, major depression and bipolar disorder. Br J Psychiatry. 2004;185:258–259. [PubMed]
  • Csernansky JG, Joshi S, Wang L, Haller JW, Gado M, Miller JP, Grenander U, Miller MI. Hippocampal morphometry in schizophrenia by high-dimensional brain mapping. Proc Natl Acad Sci USA. 1998;95:11406–11411. [PubMed]
  • Csernansky JG, Wang L, Joshi SC, Ratnanather JT, Miller MI. Computational anatomy and neuropsychiatric disease: Probabilistic assessment of variation and statistical inference of group difference, hemispheric asymmetry, and time-dependent change. Neuroimage. 2004;23:S139–S150. [PubMed]
  • Dupuis P, Grenander U, Miller MI. Variational problems on flows of diffeomorphisms for image matching. Quat Appl Math. 1998;56:587–600.
  • Fischl B, Sereno MI, Tootell RB, Dale AM. High-resolution intersubject averaging and a coordinate system for the cortical surface. Hum Brain Mapp. 1999;8:272–284. [PubMed]
  • Genovese CR, Lazar NA, Nichols T. Thresholding of statistical maps in functional neuroimaging using the false discovery rate. Neuroimage. 2002;15:870–878. [PubMed]
  • Golub G, van der Vorst HA. Eigenvalue computation in the 20th century. J Comput Appl Math. 2000;123:35–65.
  • Han X, Xu C, Prince JL. A topology preserving deformable model using level set in CVPR’2001. Vol. 2. IEEE; Kauai, HI: 2001. pp. 765–770.
  • Han X, Xu C, Braga-Neto U, Prince J. Topology correction in brain cortex segmentation using a multiscale, graph-based algorithm. IEEE Trans Med Imag. 2002;21:109–121. [PubMed]
  • Harasty J, Seldon HL, Chan P, Halliday G, Harding A. The left human speech processing cortex is thinner but longer than the right. Laterality. 2003;8:247–260. [PubMed]
  • Hirayasu Y, McCarley RW, Salisbury DF, Tanaka S, Kwon JS, Frumin M, Snyderman D, Yurgelun-Todd D, Kikinis R, Jolesz FA, Shenton ME. Planum temporale and Heschl gyrus volume reduction in schizophrenia: A magnetic resonance imaging study of first-episode patients. Arch Gen Psychiatry. 2000;57:692–699. [PMC free article] [PubMed]
  • Honeycutt NA, Musick A, Barta PE, Pearlson GD. Measurement of the planum temporale (PT) on magnetic resonance imaging scans: Temporal PT alone and with parietal extension. Psychiatry Res. 2000;98:103–116. [PubMed]
  • Joshi S, Miller MI. Landmark matching via large deformation diffeomorphisms. IEEE Trans Imag Process. 2000;9:1357–1370. [PubMed]
  • Joshi S, Miller MI, Grenander U. On the geometry and shape of brain submanifolds. Int J Pattern Recogn Artif Intell. 1997;11:1317–1343.
  • Joshi M, Cui J, Doolittle K, Joshi S, Van Essen D, Wang L, Miller MI. Brain segmentation and the generation of cortical surfaces. Neuroimage. 1999;9:461–476. [PubMed]
  • Kasai K, Shenton ME, Salisbury DF, Hirayasu Y, Onitsuka T, Spencer MH, Yurgelun-Todd DA, Kikinis R, Jolesz FA, McCarley RW. Progressive decrease of left heschl gyrus and planum temporale gray matter volume in first-episode schizophrenia: A longitudinal magnetic resonance imaging study. Arch Gen Psychiatry. 2003;60:766–775. [PMC free article] [PubMed]
  • Kessler RC, Andrews G, Mroczek D, Ustun B, Wittchen H-U. The world health organization composite international diagnostic interview short form (CIDI-SF) Int J Methods Psychiatr Res. 1998;7:33–55.
  • Kwon J, McCarley RW, Hirayasu Y, Anderson J, Fischer IA, Kikinis R, Jolesz FA, Shenton M. Left planum temporale volume reduction in schizophrenia. Arch Gen Psychiatry. 1999;56:142–148. [PubMed]
  • Miller MI, Massie AB, Ratnanather JT, Botteron KN, Csernansky JG. Bayesian construction of geometrically based cortical thickness metrics. Neuroimage. 2000;12:676–687. [PubMed]
  • Miller MI, Hosakere M, Barker AR, Priebe CE, Lee N, Ratnanather JT, Wang L, Gado M, Morris JC, Csernansky JG. Labeled cortical mantle distance maps of the cingulate quantify differences between dementia of the Alzheimer type and healthy aging. Proc Natl Acad Sci USA. 2003;100:15172–15177. [PubMed]
  • Miller MI, Beg MF, Ceritoglu C, Stark C. Increasing the power of functional maps of the medial temporal lobe by using large deformation diffeomorphic metric mapping. Proc Natl Acad Sci USA. 2005;102:9685–9690. [PubMed]
  • Nakagawa A. Role of anterior and posterior attention networks in hemispheric asymmetries during lexical decisions. J Cogn Neurosci. 1991;3:315–321. [PubMed]
  • Nieto-Castanon A, Ghosh SS, Tourville JA, Guenther FH. Region of interest based analysis of functional imaging data. Neuroimage. 2003;19:1303–1316. [PubMed]
  • Pentland A, Sclaroff S. Closed-form solutions for physically based shape modeling and recognition. IEEE Trans Pattern Anal Mach Intell. 1991;12:715–729.
  • Qiu A. PhD Thesis. 2006. Intrinsic and Extrinsic Analysis in Computational Anatomy. [PubMed]
  • Qiu A, Bitouk D, Miller MI. Smooth functional and structural maps on the neocortex via orthonormal bases of the Laplace–Beltrami operator. IEEE Trans Med Imag. 2006;25:1296–1306. [PubMed]
  • Ratnanather JT, Barta PE, Honneycutt NA, Lee N, Morris NG, Dziorny AC, Hurdal MK, Pearlson GD, Miller MI. Dynamic programming generation of boundaries of local coordinatized submanifolds in the neocortex: Application to the planum temporale. Neuroimage. 2003;20:359–377. [PubMed]
  • Robb RA, Hanson DP, Karwoski RA, Larson AG, Workman EL, Stacy MC. Analyze: A comprehensive, operator-interactive software package for multidimensional medical image display and analysis. Comput Med Imag Graph. 1989;13:433–454. [PubMed]
  • Seldon HL. Does brain white matter growth expand the cortex like a balloon? Hypothesis and consequences. Laterality. 2005;10:81–95. [PubMed]
  • Shapleske J, Rossell SL, Simmons A, David AS, Woodruff PW. Are auditory hallucinations the consequence of abnormal cerebral lateralization? A morphometric MRI study of the sylvian fissure and planum temporale. Biol Psychiatry. 2001;49:685–693. [PubMed]
  • Sheehan DV, Lecrubier Y, Sheehan KH, Amorim P, Janavs J, Weiller E, Hergueta T, Baker R, Dunbar GC. The Mini-International Neuropsychiatric Interview (M.I.N.I.): The development and validation of a structured diagnostic psychiatric interview for DSM-IV and ICD-10. J Clin Psychiatry. 1998;59:22–23. [PubMed]
  • Thompson PM, Hayashi KM, Sowell ER, Gogtay N, Giedd JN, Rapoport JL, de Zubicaray GI, Janke AL, Rose SE, Semple J, Doddrell DM, Wang Y, van Erp TG, Cannon TD, Toga AW. Mapping cortical change in Alzheimer’s disease, brain development, and schizophrenia. Neuroimage. 2004;23:S2–S18. [PubMed]
  • Tosun D, Rettmann ME, Han X, Tao X, Xu C, Resnick SM, Pham DL, Prince JL. Cortical surface segmentation and mapping. Neuroimage. 2004a;23:S108–S118. [PubMed]
  • Tosun D, Rettmann ME, Prince JL. Mapping techniques for aligning sulci across multiple brains. Med Image Anal. 2004b;8:295–309. [PubMed]
  • Trouvé A. An infinite dimensional group approach for physics based models. Technical report. 1995. (electronically available at Unpublished.
  • Vaillant M, Glaunès J. Surface matching via currents. Lecture Notes Comp Sci: Inform Proc Med Imag. 2005;3565:381–392. [PubMed]
  • Vaillant M, Qiu A, Glaunès J, Miller MI. Diffeomorphic metric surface mapping in subregion of the superior temporal gyrus. Neuroimage. 2007;34:1149–1159. [PMC free article] [PubMed]
  • Van Essen D. Surface-based approaches to spatial localization and registration in primate cerebral cortex. Neuroimage. 2004;23:S97–S107. [PubMed]
  • Van Essen D. A population-average, landmark- and surface-based (pals) atlas of human cerebral cortex. Neuroimage. 2005;28:635–662. [PubMed]
  • Van Essen D, Drury H. Structural and functional analyses of human cerebral cortex using a surface-based atlas. J Neurosci. 1997;17:7079–7102. [PubMed]
  • Van Essen DC, Drury HA, Joshi S, Miller MI. Functional and structural mapping of human cerebral cortex: Solutions are in the surfaces. Proc Natl Acad Sci USA. 1998;95:788–795. [PubMed]
  • Van Essen D, Drury H, Dickson J, Harwell J, Hanlon D, Anderson C. An integrated software suite for surface-based analyses of cerebral cortex. J Am Med Inform Assoc. 2001;8:443–459. [PMC free article] [PubMed]
  • von Economo C. The Cytoarchitectonics of the Human Cerebral Cortex. Oxford: Oxford University Press; 1929.
  • Weisbrod M, Maier S, Harig S, Himmelsbach U, Spitzer M. Lateralized semantic and indirect semantic priming effects in people with schizophrenia. Br J Psychiatry. 1998;172:142–146. [PubMed]
  • Worsley K. Local maxima and the expected Euler characteristic of excursion sets of χ2, F and t fields. Adv Appl Probab. 1994;26:13–42.
  • Worsley KJ, Marrett S, Neelin P, Vandal AC, Friston KJ, Evans AC. A unified statistical approach for determining significant signals in images of cerebral activation. Hum Brain Mapp. 1996;4:58–73. [PubMed]
  • Worsley K, Taylor J, Tomaiuolo F, Lerch J. Unified univariate and multivariate random field theory. Neuroimage. 2004;23:S189–S195. [PubMed]
  • Zatorre RJ, Belin P, Penhune VB. Structure and function of auditory cortex: Music and speech. Trends Cogn Sci. 2002;6:37–46. [PubMed]