Detection of fMRI activation using Cortical Surface Mapping ^†

General Overview

CSM methodology can be subdivided into several stages. The procedures involved, which are briefly described below, are either common to standard, 3D‐oriented methods or require adaptations, to a greater or lesser extent. Important conceptual differences with respect to standard methods are presented in a separate treatment of the procedures in question. The stages comprise:

• Spatial image transformation(s): The set of functional MRI scans acquired in the course of the experimental protocol can be spatially realigned and/or normalised to match a template. These are optional spatial transformations, and the decision to apply one or both should be guided by criteria (subject motion, purpose of making inter‐subject comparisons) that have no direct bearing upon the specific nature of the proposed methodology. Coregistration between the T${}_{\text{*}}^{\text{2}}$‐weighted scans and a T₁‐weighted anatomical scan can be performed to ensure a satisfactory match between the cortical representation (obtained from the T₁‐weighted image) and the functional volumes. The assignment of values to cortical positions (see below) relies heavily upon this matching; therefore this is a critical step. However, coregistring is not compulsory, provided that no substantial subject movement takes place between the functional and anatomical acquisitions. Distortion in the functional scans, due to Echo‐Planar acquisition, is also a major issue, since it can put in jeopardy the possibility of achieving a satisfactory match. For moderate field MRI machines, distortion should not be strong enough to require the application of a correction algorithm [Jezzard and Clare, 1999].

• Grey matter/white matter segmentation: Grey and white matter masks are built through segmentation of the T₁‐weighted anatomical scan of the subject [Mangin et al., 1995].

• Surface extraction: A triangulated lattice representation of the internal surface of the cortex is extracted, through application of a surface tracking algorithm to the white matter segmented volume [Gordon and Udupa, 1989]. The lattice comes under the form of a very dense mesh (the node separation is of the order of the anatomical acquisition voxel size, typically ∼1 mm), with a size that varies typically between 100,000 and 150,000 nodes for each hemisphere. Constraints in the segmentation method insure that the lattice is homotopic to a sphere.

• Value assignment to the surface positions: The T${}_{\text{*}}^{\text{2}}$‐weighted functional scans are interpolated on the positions corresponding to the nodes of the cortical lattice (plus a user defined outward shift in the surface normal direction to account for cortical thickness). The outcome of this operation is a set of cortical surface scans (CS‐scans) with n values, n standing for the number of nodes in the cortical lattice. Each value is associated with a set of coordinates corresponding to the interpolating position. Subsequent analysis steps will thus be confined to a set of values interpolated from the neighbourhood of the cortex.

• Spatial smoothing: Each one of the CS‐scans, consisting of a set of values associated with node positions, is spatially smoothed. While the purpose behind this surface‐based smoothing is essentially the same as for the classical 3D analysis, the conceptual basis of its implementation is markedly different, and therefore this step is treated in more detail below.

• General Linear Model parameter estimation: Estimation is performed in essentially the same way as with classical three‐dimensional analysis. The option to perform global correction (proportional scaling) or temporal filtering is retained. The output is a set of estimated parameters, with the same dimension (i.e., number of nodes) as an individual functional CS‐scan.

• Smoothness estimation: For this step, we have employed “statistical flattening” [Worsley et al., 1999], a generalised estimation procedure that takes into account local fluctuations of smoothness, and is easily adaptable to the particular requirements of the proposed methodology. This procedure and its implementation receive a more thorough (albeit informal) independent treatment below.

• Hypothesis testing based on intensity: The assignment, for each node, of a P value reflecting the risk of error concerning the rejection of a null hypothesis involving a linear combination of the estimated parameters is based upon the theory of random fields, and does not differ from currently employed procedures [Friston et al., 1995; Worsley et al., 1996b].

Surface‐based Smoothing

The purpose behind spatial smoothing of functional brain imaging data is two‐fold [Petersson et al., 1999]:

• Improving signal‐to‐noise ratio, through the use of a smoothing kernel wide enough to filter out a significant part of the (mostly high‐frequency) noise.

• Ensuring that the output field of values constitutes an adequate discrete representation of a continuous statistical field, so that subsequent inference procedures based on the theory of random fields are valid. For this to be achieved, the smoothness of the field must be large compared to the sampling resolution of the images under study. Note that this is not an absolute constraint, since subsampling can always be employed as a way to impose the required conditions.

Smoothing is usually implemented as a simple, fully isotropic convolution with a three‐dimensional Gaussian kernel, and it is applied uniformly over the entire volume. Such an implementation does not take into account the anatomical configuration of the brain (different kinds of tissue, deeply folded structure, etc.). This leads to averaging of signal emanating from the cortex and from while matter of CSF, as well as to averaging of signal issuing from sources that are close to each other in a Euclidean sense, but geodesically far apart. CSM methodology implements surface‐oriented smoothing, aiming to obtain a spatial correlation structure that depends on geodesic rather than Euclidean distance. In this way, the introduction of artificial correlation between areas that are close together in Euclidean terms but geodesically far apart (e.g., two points on opposite sides of a sulcus) is avoided.

This implementation of surface‐oriented smoothing differs from simple, voxel‐based isotropic three‐dimensional convolution in two main aspects:

• At the practical level: smoothing is applied to a field of values sampled with an irregular grid (as opposed to a regular voxel‐based sampling).

• At the conceptual level: smoothing is done along the cortical surface, instead of being extended to the whole acquired volume.

Clearly, convolution with a Gaussian kernel is no longer a viable choice: the highly convoluted nature of the cortex, along with the irregular sampling, requires a new theoretical framework to build the smoothing procedure upon.

Smoothing an irregularly sampled field of values

The solution adopted to take the irregularity of the sampling grid into account relies on a result widely employed in image processing [Koenderink, 1984]: the equivalence between Gaussian convolution and heat diffusion. The diffusion of heat in a time‐varying field of temperatures I(r→, t) proceeds according to a law expressed as

(1)

where t is time, K is a constant involving several physical parameters (K = c/ρc_p, with c standing for thermal conductivity, ρ for density and c _p for specific heat), and ∇² is the Laplacian operator (∇² I(x, y) = ∂²I/∂x² + ∂²I/∂y², in two dimensions). This assumes that K is constant all over the field. On the other hand, if t is taken to represent a smoothness (scale‐space) parameter, instead of time, it can be proved that

(convolution of an image I(r→,t = 0) with a Gaussian kernel) is a solution of Equation 1 (with K = 0.5). Hence, applying a Gaussian spatial filter with parameter to an image is the same thing as letting a diffusion process run its course for a lapse of time T. The advantage of this formulation based on a differential equation is that it lends itself more easily to being adapted for the specific case of an irregular 2D lattice embedded in a three‐dimensional space.

Numerical computation of the results of a diffusion process over a discrete field of values can be carried out as an iterative process of the form

(2)

for each temporal iteration step (Δt). is the local estimation of the Laplacian at time t. The application of this solution would be trivial in the case of a regular lattice with a fixed neighbourhood: finite differences, for instance, could be chosen as a method for locally estimating partial derivatives. For an irregular lattice, a general approach can be employed [Liszka and Orkisz, 1980; Huiskamp, 1991], based on the Taylor series expansion of a function around a point. Assume that an image, implicitly described by a function I(x, y), is sampled with an irregular lattice composed of interconnected nodes. For a point (x₀,y₀) of the image, this expansion has the form:

(3)

where I_i = I(x_i, y_i), h_i = x_i − x₀, k_i = y_i − y₀, and . Writing Equation 3 for a surface element consisting of a lattice node, located at (x ₀, y ₀), and its neighbours i = 1, 2,…m, we obtain the set of equations

where

and the five derivatives at point (x₀, y₀) are

In this fashion, estimates of the two‐dimensional Laplacian (∇² I(x, y) = ∂² I/∂x² + ∂² I/∂y²) are obtained at each node of the lattice. This approach amounts to solving, for each node, a linear system involving the relative positions of the node and its neighbours and the corresponding field value differences. It is still a finite differences method, in the sense that the partial derivatives are estimated by differences between field values in neighbouring points, but its range of application is extended to any arbitrarily irregular grid. Some constraints remain, however, concerning the configuration of the lattice. Firstly, the application of post‐processing procedures (e.g., decimation) to the lattice may lead to situations in which node neighbours are fewer than the number of variables estimated locally (this number will be 5 in the case of the second order Taylor expansion necessary to estimate the Laplacian, or 4 if isotropy is assumed and the ∂²f₀/∂x∂y term is considered to be zero), making the corresponding linear system underdetermined, which means that the solution will not be unique. One of the options to pick one among the infinity of solutions thus obtained is to select the solution with the smallest norm. This is the condition implied by the pseudo‐inverse based resolution that we adopted. This was shown to be a sensible option in practice: provided that the nodes with less than four neighbors constitute a small (∼1%) proportion of the total number of nodes, the final results do not differ significantly from those obtained in a situation of system determinacy all over the lattice. The second constraint has to do with node separation: the nodes must be sufficiently close together for the partial derivative estimations to be reliable. The specifics of lattice building are not discussed in this paper, but it should be made clear that these constraints can be applied at the surface extraction stage, by suitable tuning of the extraction algorithm parameters. In the analyses we carried out, we employed both undecimated lattices (the direct output of the surface extraction stage), with ∼100,000–150,000 nodes per hemisphere, and lattices that underwent constrained node decimation. The latter had ∼25,000 nodes per hemisphere, and a small (<1%) proportion of nodes with less than four neighbours. Mean and maximal internodal distance were, respectively, 1 and 1.8 mm (undecimated lattice) and 2.5 and 6 mm (decimated lattice). Only minor differences were detected between the results obtained with undecimated and decimated lattices. This seems to indicate that the presence of important internodal distance variability, and the resulting local decrease in the accuracy of Laplacian estimation, does not per se affect the outcome of the analysis in a significant way.

Another important concern is the choice of the temporal iteration step (Δt in Equation 2). The critical value above which the numerical resolution will be unstable [Press et al., 1992; Gerig et al., 1992], in the specific case of the diffusion equation, depends on thermal conductivity, a parameter included in the K term in Equation 1. In the case of diffusion‐based smoothing over a lattice, this conductivity term can be seen as conceptually equivalent to node distance: the closer a pair of nodes is, the stronger their mutual influence, and this mimics the situation of high thermal conductivity and fast heat propagation. Irregular node spacing, thus, parallels a situation of varying thermal conductivity across the field. The critical temporal iteration step is directly proportional to the minimal node distance, because instability tends to spread over the whole lattice, and therefore a single pair of nodes that are too close together will be enough to globally affect the results. Lattices possessing nodes in these conditions require very small temporal iteration steps, and this may lead to exceedingly long processing times. Hence, a double constraint operates at the level of node separation: small internodal distances imply a small temporal iteration step and long computation times; on the other hand, very large internodal distances will eventually affect the reliability of the Laplacian estimations. There is, thus, a trade‐off between processing speed and estimation accuracy.

If computation time is not unreasonably large, a conservative value for Δt should be chosen. We have adopted a default value of Δt = 0.1, throughout the analysis of the fMRI protocol described below, which proved to be low enough to lead to a stable process for the undecimated lattices we worked with (minimal internodal distance ∼0.6–0.7 mm). For these very dense lattices (∼100,000–150,000 nodes), this choice implies that a diffusion process equivalent to convolving with a Gaussian kernel with FWHM = 8 mm (σ = 3.397) requires 115 iterations (that is σ²/Δt rounded to the nearest integer). In terms of computation time, this represents roughly 2–3 min in a SUN ULTRA 10 workstation, for a single CS‐scan (one hemisphere). Grouping scans in blocks treated simultaneously allows to optimise smoothing time, so that the increase in time is not linear with respect to the total number of scans. 8 mm smoothing for a protocol consisting of ∼100 scans takes a few hours, making this by far the most time‐consuming step of the analysis. Computation time may quickly rise to prohibitive values for bulkier protocols, and this is why the use of decimated lattices may be recommended, or even become compulsory. To give an idea of the decrease in computation time, the same amount of smoothing applied to the same protocol, but this time with scans built from decimated lattices (∼25,000 nodes) took ∼40 min.

Smoothing Along a Curved Surface

The theoretical framework described in the previous paragraph was chosen to tackle the problem of smoothing an irregularly sampled field of values. We have focused on its 2D version, since the aim is to smooth over the surface. However, the fact that the cortical representation lattice is embedded in a volume raises a problem concerning the proper axis system upon which to base the estimation of the local partial derivatives. A conflict exists between the need to employ three spatial coordinates in the Euclidean space to unambiguously describe the position of each node and the wish to perform a strictly surface‐based smoothing: locally estimating three‐dimensional Laplacians and solving the diffusion equation accordingly would clearly go against this intent. The definition of a coordinate system intrinsic to the surface would allow for a strictly 2D‐based smoothing. This implies a parameterization, i.e., the definition of a mapping function m such that (x, y, z) ≡ m⁻¹ (u, v), u and v being the new coordinates that allow to refer to every point in the surface without ambiguity. The diffusion equation is then solved along u and v, in a strictly 2D fashion.

The parameterization adopted was a simple local transformation that maps each surface element (a node and its first neighbours) into the plane, while keeping unchanged both the edge distances and the angular proportions between the edges [Welch and Witkin, 1994]. An individual mapping function m _i is defined for each node i, independently of the surrounding surface elements. This approach is made possible since, for each iteration step of the numerical solution of the diffusion equation, the estimations involve only the differences between the values associated with each node and its nearest neighbours. Locally mapping each surface element into the plane avoids the severe areal distortion that would result from a global flattening. This distortion is of the order of 10–20% [Drury et al., 1996], but locally attains much higher values.

Multiple Comparisons Correction by Means of “Statistical Flattening”

Statistical inference of parametric maps must take into account the degree of effective mutual dependence presented by the individual hypothesis tests. Considering the tests as totally independent would not be realistic, since, due both to physiological reasons and to instrumentation‐related specifics of fMRI acquisition, a certain amount of correlation will always be present, and hence a Bonferroni‐type correction would generally prove too harsh [Worsley et al., 1992]. Usually, an estimation of the field smoothness is made, and the stringency of the correction depends upon this measure of spatial correlation.

Smoothness‐based correction can be seen as the computation of the number of tests normalised for the global smoothness of the statistical field. This can alternatively be seen as “shrinking” the field so that its extent, in each of the space directions, is expressed in smoothness‐normalised Resolution Elements (Resels) [Worsley et al., 1992] instead of voxels. An intuitive generalisation of this principle consists in the estimation of local smoothness parameters, and in locally shrinking the field so that distances are normalised for the smoothness estimation at that point. This, in very general terms, is the idea behind the “statistical flattening” correction procedure, described elsewhere [Worsley et al., 1999], and chosen to implement multiple comparisons correction in the cortical surface mapping methodology described herein. This procedure applies local smoothness estimations based on the residual values, after removal of the effects accounted for by the experimental linear model. To permit exact flattening, this operation takes place in the multidimensional residual space, impossible to visualize, into which the individual values are mapped.

Although statistical flattening can be applied in the case of standard volumic analyses (it has, indeed, been incorporated to the latest release of the SPM package), it lends itself particularly well to being adapted for the case of data associated with a triangulated lattice: edge distances are normalised by the correlation between the normalised residuals associated with the corresponding nodes, so that the ensuing correlation structure is isotropic, and the total surface of the deformed field reflects the number of effectively independent tests involved.