Analysis of fMRI data using cortical surface models offers at least three advantages over more conventional 3-dimensional analysis methods. First, cortical surface models allow for better visualization of activations, providing a more global view than single slices and a better view of the spatial extent of activation foci and their locations relative to each other and to sulcal/gyral landmarks (

Dale and Sereno, 1993). Second, statistical methods for the analysis of single subject data can benefit from the exclusion of non-gray matter signals, and smoothing signals along the cortical surface, rather than in 3D results in superior resolution and sensitivity (

Kiebel et al., 2000;

Andrade et al., 2001;

Formisano et al., 2004). Finally, group analysis with cortical surface models employs inter-subject alignment based on the patterns of sulci and gyri, as opposed to Talairach registration which often ignores sulcal/gyral landmarks and tends to blur activity across neighboring banks of a sulcus (

Fischl et al., 1999a;

Fischl et al., 1999b).

In this paper we describe methods that we have used and developed to facilitate cortical surface-based inter-subject analyses. A routine aspect of inter-subject fMRI analyses -- for both 2D cortical surface and 3D volumetric methods -- is spatial smoothing. Smoothing acts as a low-pass spatial frequency filter and thus improves the signal-to-noise ratio (SNR) by filtering out high spatial frequency noise (

Petersson et al., 1999a). Smoothing increases the likelihood that inter-subject analyses will detect signals from foci that display variability in their precise cortical location across subjects. Spatial smoothing also ensures that the data more closely approximate a continuous field of random values, a necessary assumption of the random field theory used for multiple comparisons correction (

Worsley et al., 1992;

Worsley et al., 1996;

Petersson et al., 1999b;

Worsley et al., 1999;

Andrade et al., 2001).

The vertices of a cortical surface mesh are not, however, arranged on a grid with regular spacing; instead distances and angles between neighboring vertices are slightly variable (

Fischl et al., 1999a;

Chung et al., 2003). Because of this, direct application of a Gaussian blurring kernel -- as is done with 3D datasets -- is computationally intensive for surface-based data. A computationally efficient method is to iteratively perform nearest-neighbor averaging; where each vertex's value is averaged with those of its neighbors. Slightly more complicated iterative smoothing algorithms have also been developed, which rely on a heat diffusion model and involve unequally weighting the contribution from each neighbor in the post-iteration value of a given vertex (

Andrade et al., 2001;

Chung et al., 2003;

Chung et al., 2005). For example, in the latest and simplest version of this method, called heat kernel smoothing, the weights are calculated based on the distances between a vertex and each of its neighbors (

Chung et al., 2005). The effect of such iterative smoothing is quite similar to what would be expected if a Gaussian blurring kernel were used.

For heat kernel smoothing, the full-width-half-max (FWHM) size of the Gaussian blurring kernel equivalent to a particular number of iterations can be predicted from heat diffusion equations (

Chung et al., 2003;

Chung et al., 2005). The accuracy of those predictions, however, is limited by the extent to which the assumptions of the heat diffusion model are met. Specifically, the model assumes continuous diffusion, whereas iterative smoothing is inherently discrete and the vertices are spaced at discrete intervals. The bandwidth of the heat kernel determines the amount of smoothing performed at each iteration, and larger bandwidths make the continuous diffusion model less accurate (

Chung et al., 2005). Similarly, the heat diffusion model is less accurate with increasing inter-neighbor distances. Thus, even though the heat diffusion model provides predictions of FWHM smoothness as a function of number of iterations and bandwidth, it is necessary to determine the accuracy of those predictions. To address this concern, we have made empirical measures of FWHM smoothness as a function of smoothing steps for both heat kernel smoothing and nearest-neighbor average smoothing.

Another important aspect of inter-subject analyses is the statistical correction for multiple comparisons. Multiple comparisons correction is desirable because of the large number of voxels, and thus independent statistical tests, involved in voxel-by-voxel analysis of fMRI images. An alternative to the overly conservative Bonferroni correction is cluster size exclusion (

Worsley et al., 1992;

Poline and Mazoyer, 1993;

Friston et al., 1994;

Forman et al., 1995;

Worsley et al., 1999;

Andrade et al., 2001). Modeling the value (e.g. t-statistic) of each voxel or vertex as a normally distributed random variable, it is assumed that neighbors display some degree of covariance, or spatial smoothness. That is, the value at a voxel is likely to be somewhat similar to its neighbors, either because of spatial smoothing, spatially correlated noise -- produced by the scanner or physiologically -- or because of the spatial extent of actual brain activations and the resulting hemodynamic response (

Kiebel et al., 1999). With smoother data, it is more likely that a suprathreshold voxel will have neighbors that are also suprathreshold, forming a cluster.

Cluster size limits can be generated, with random field theory (

Worsley et al., 1996;

Andrade et al., 2001), permutation tests (

Hayasaka and Nichols, 2003,

2004;

Hayasaka et al., 2004), or Monte Carlo simulations (

Poline and Mazoyer, 1993;

Forman et al., 1995), for any number of uncorrected p-values. A liberal p-value coupled with a very large cluster size limit is theoretically equivalent to a conservative p-value coupled with a small cluster size limit. In practice, however, the choice of the uncorrected p-value can have a strong influence on the number of clusters satisfying the cluster size limits. Furthermore, for the purpose of defining regions of interest (ROIs), the use of a conservative threshold will tend to underestimate the true size of a focus of activation, identifying only the very peak of activations; lower thresholds will tend to result in the joining of distinct clusters. This situation has motivated our development of a method for identifying ROIs using a sliding threshold followed by growth of clusters.

Because random cluster sizes depend on the smoothness of the images, an estimate of this is necessary. Actual data, however, contains regions of activation with some spatial extent, increasing measures of overall smoothness which would bias the cluster size thresholds to be overly conservative (

Kiebel et al., 1999). For this reason, it would be more appropriate to measure the smoothness of random noise that shares the spatial correlations related to hemodynamic and scanner properties. As a proxy for actually measuring noise, smoothness can be measured from the normalized residual error of single subject GLM deconvolution (

Kiebel et al., 1999). In order to determine a suitable method for estimating the intrinsic smoothness of a group analysis dataset, we compared the smoothness of actual group averages, normalized residual error of group analysis, and the group average of simulated noise datasets.