Statistical Analysis on Brain Surfaces

A local parameterization is usually done by fitting a quadratic polynomial of the form

$\displaystyle{\bf X}(u^{1},u^{2})=\beta_{1}u^{1}+\beta_{2}u^{2}+\beta_{3}(u^{1})^{2}+\beta_{4}u^{1}u^{2}+\beta_{5}(u^{2})^{2}$

(1)

in $(u^{1},u^{2})\in D\subset\mathbb{R}^{2}$. The data can be centered so there is no constant term in the quadratic form (1) (Chung, 2012). The coefficients $\beta_{i}$ are usually estimated by the least squares method (Joshi et al., 1995; Khaneja et al., 1998; Chung, 2012).

In estimating various differential geometric measures such as the Laplace-Beltrami operator or curvatures, it is not necessary to find global surface parameterization of $\mathcal{M}$. Local surface parameterization such as the quadratic polynomial fit is sufficient to obtain such geometric quantities (Chung, 2012). The drawback of the polynomial parameterization is that there is a tendency to weave the outermost mesh vertices to find vertices in the center. Therefore this is not advisable to directly fit (1) when one of the coordinate values rapidly changes.

Parameterizing cortical and subcortical surfaces with respect to simpler algebraic surfaces such as a unit sphere is needed to establish a standard coordinate system. However, polynomial regression type of local parameterization is ill-suited for this purpose. For the global surface parameterization, we can use surface flattening (Andrade et al., 2001; Angenent et al., 1999), which is nonparametric in nature. The surface flattening parametrizes a surface by either solving a partial differential equation or optimizing its variational form.

Deformable surface algorithms naturally provide one-to-one maps from cortical surfaces to a sphere since the algorithm initially starts with a spherical mesh and deforms it to match the tissue boundaries (MacDonald et al., 2000). The deformable surface algorithms usually start with the second level of triangular subdivision of an icosahedron as the initial surface. After several iterations of deformation and triangular subdivision, the resulting cortical surface contains very dense triangle elements. There are many surface flattening techniques such as conformal mapping (Angenent et al., 1999; Gu et al., 2004; Hurdal & Stephenson, 2004) quasi-isometric mapping (Timsari & Leahy, 2000), area preserving mapping (Brechbuhler et al., 1995), and the Laplace equation method (Chung, 2012).

For many surface flattening methods to work, the starting binary object has to be close to star-shape or convex. For shapes with a more complex structure, the methods may create numerical singularities in mapping to the sphere. Surface flattening can destroy the inherent geometrical structure of the cortical surface due to the metric distortion. Any structural or functional analysis associated with the cortex can be performed without surface flattening if the intrinsic geometric method is used (Chung, 2012).

The spherical harmonic (SPHARM) representation ²²2The SPHARM package is available at
www.stat.wisc.edu/$\sim$mchung/softwares/weighted-SPHARM/weighted-SPHARM.html is a widely used subcortical surface parameterization technique (Chung et al., 2008; Chung, 2012; Gerig et al., 2001; Gu et al., 2004; Kelemen et al., 1999; Shen et al., 2004). SPHARM represents the coordinates of mesh vertices as a linear combination of spherical harmonics. SPHARM has been mainly used as a data reduction technique for compressing global shape features into a small number of coefficients. The main global geometric features are encoded in low degree coefficients while the noise are in high degree spherical harmonics (Gu et al., 2004). The method has been used to model various subcortical structures such as ventricles (Gerig et al., 2001), hippocampi (Shen et al., 2004) and cortical surfaces (Chung et al., 2008). The spherical harmonics have a global support. So the resulting spherical harmonic coefficients contain the global shape features and it is not possible to directly obtain local shape information from the coefficients only. However, it is still possible to obtain local shape information by evaluating the representation at each fixed vertex, which gives the smoothed version of the coordinates of surfaces. In this fashion, SPHARM can be viewed as mesh smoothing (Chung et al., 2008, n.d.). In this section, we present a brief introduction of SPHARM within a Hilbert space framework.

Suppose there is a bijective mapping between the cortical surface $\mathcal{M}$ and a unit sphere $S^{2}$ obtained through a deformable surface algorithm (Chung, 2012). Consider the parameterization of $S^{2}$ by

$${\bf X}(\theta,\varphi)=(\sin\theta\cos\varphi,\sin\theta\sin\varphi,\cos\theta),$$

with $(\theta,\varphi)\in[0,\pi)\otimes[0,2\pi)$. The polar angle $\theta$ is the angle from the north pole and the azimuthal angle $\varphi$ is the angle along the horizontal cross-section. Using the bijective mapping, we can parameterize functional data $f$ with respect to the spherical coordinates

$\displaystyle f(\theta,\varphi)=g(\theta,\varphi)+\epsilon(\theta,\varphi),$

(2)

where $g$ is a unknown smooth coordinate function and $\epsilon$ is a zero mean random field, possibly Gaussian. The error function $\epsilon$ accounts for possible mapping errors. The unknown signal $g$ is then estimated in the finite subspace of $\mathcal{L}^{2}(S^{2})$, the space of square integrable functions in $S^{2}$, spanned by spherical harmonics in the least squares fashion (Chung et al., 2008).

Previous imaging and shape modeling literature have used the complex-valued spherical harmonics (Bulow, 2004; Gerig et al., 2001; Gu et al., 2004; Shen et al., 2004). In practice, however, it is sufficient to use only real-valued spherical harmonics (Courant & Hilbert, 1953; Homeier & Steinborn, 1996), which is more convenient in setting up a real-valued stochastic model (2). The relationship between the real- and complex-valued spherical harmonics is given in (Blanco et al., 1997; Homeier & Steinborn, 1996). The complex-valued spherical harmonics can be transformed into real-valued spherical harmonics using an unitary transform.

The spherical harmonic $Y_{lm}$ of degree $l$ and order $m$ is defined as

$\displaystyle Y_{lm}=\begin{cases}c_{lm}P^{|m|}_{l}(\cos\theta)\sin(|m|\varphi),&-l\leq m\leq-1,\\ \frac{c_{lm}}{\sqrt{2}}P_{l}^{|m|}(\cos\theta),&m=0,\\ c_{lm}P^{|m|}_{l}(\cos\theta)\cos(|m|\varphi),&1\leq m\leq l,\end{cases}$

where $c_{lm}=\sqrt{\frac{2l+1}{2\pi}\frac{(l-|m|)!}{(l+|m|)!}}$ and $P^{m}_{l}$ is the associated Legendre polynomial of order $m$ (Courant & Hilbert, 1953; Wahba, 1990), which is given by

$\displaystyle P_{l}^{m}(x)=\frac{(1-x^{2})^{m/2}}{2^{l}l!}\frac{d^{l+m}}{dx^{l+m}}(x^{2}-1)^{l},x\in[-1,1].$

The first few terms of the spherical harmonics are

	$\displaystyle Y_{00}=\frac{1}{\sqrt{4\pi}},Y_{1,-1}=\sqrt{\frac{3}{4\pi}}\sin\theta\sin\varphi,$
	$\displaystyle Y_{1,0}=\sqrt{\frac{3}{4\pi}}\cos\theta,Y_{1,1}=\sqrt{\frac{3}{4\pi}}\sin\theta\cos\varphi.$

The spherical harmonics are orthonormal with respect to the inner product

$\displaystyle\langle f_{1},f_{2}\rangle$

$\displaystyle=$

$\displaystyle\int_{S^{2}}f_{1}(\Omega)f_{2}(\Omega)\;d\mu(\Omega),$

where $\Omega=(\theta,\varphi)$ and the Lebesgue measure $d\mu(\Omega)=\sin\theta d\theta d\varphi$. The norm is then defined as

$\displaystyle||f_{1}||=\langle f_{1},f_{1}\rangle^{1/2}.$

(3)

The unknown mean function $g$ is estimated by minimizing the integral of the squared residual in $\mathcal{H}_{k}$, the space spanned by up to $k$-th degree spherical harmonics:

$\displaystyle\widehat{g}(\Omega)=\arg\min_{h\in\mathcal{H}_{k}}\int_{S^{2}}\Big{|}f(\Omega)-h(\Omega)\Big{|}^{2}\;d\mu(\Omega).$

(4)

It can be shown that the minimization is given by

$\displaystyle\widehat{g}(\Omega)=\sum_{l=0}^{k}\sum_{m=-l}^{l}\langle f,Y_{lm}\rangle Y_{lm}(\Omega),$

(5)

the Fourier series expansion. The expansion (5) has been referred to as the spherical harmonic representation (Chung et al., 2008; Gerig et al., 2001; Gu et al., 2004; Shen et al., 2004; Shen & Chung, 2006). This technique has been used in representing various brain subcortical structures such as hippocampi Shen et al. (2004) and ventricles (Gerig et al., 2001) as well as the whole brain cortical surfaces (Chung et al., 2008; Gu et al., 2004). By taking each component of Cartesian coordinates of mesh vertices as the functional signal $f$, surface meshes can be parameterized as a function of $\theta$ and $\varphi$.

The spherical harmonic coefficients can be estimated in least squares fashion. However, for an extremely large number of vertices and expansions, the least squares method may be difficult to directly invert large matrices. Instead, the iterative residual fitting (IRF) algorithm (Chung et al., 2008) can be used to iteratively estimate the coefficients by partitioning the larger problem into smaller subproblems. The IRF algorithm is similar to the matching pursuit method (Mallat & Zhang, 1993). The IRF algorithm was developed to avoid the computational burden of inverting a large linear problem while the matching pursuit method was originally developed to compactly decompose a time-frequency signal into a linear combination of a pre-selected pool of basis functions. Although increasing the degree of the representation increases the goodness-of-fit, it also increases the number of estimated coefficients quadratically. So it is necessary to stop the iteration at the specific degree $k$, where the goodness-of-fit and the number of coefficients balance out. This idea was used in determining the optimal degree of SPHARM (Chung et al., 2008).

The limitation of SPHARM is that it produces the Gibbs phenomenon, i.e., ringing artifacts, for discontinuous and rapidly changing continuous measurements (Chung et al., 2008; Gelb, 1997). The Gibbs phenomenon can be effectively removed by weighting the spherical harmonic coefficients exponentially smaller, which makes the representation smooth out rapidly changing signals. The weighted version of SPHARM is related to isotropic diffusion smoothing (Andrade et al., 2001; Cachia, Mangin, Riviére, Kherif, Boddaert, Andrade, Papadopoulos-Orfanos, Poline, Bloch, Zilbovicius, Sonigo, Brunelle & Régis, 2003; Chung, 2012; Chung et al., 2005) as well as the diffusion wavelet transform (Chung et al., 2008, n.d.; Hosseinbor et al., 2014).

Anatomical objects extracted from 3D medical images are aligned using affine transformations to remove the global size differences. Affine registration requires identifying corresponding landmarks either manually or automatically. The affine transform $T$ of point $p=(p_{1},\cdots,p_{d})^{\prime}\in\mathbb{R}^{d}$ to $q=(q_{1},\cdots,q_{d})^{\prime}$ is given by

$$q=Rp+c,$$

where the matrix $R$ corresponds to rotation, scaling and shear while $c$ corresponds to translation. Although the affine transform is not linear, it can be made into a linear form by augmenting the transform. The affine transform can be rewritten as

$\displaystyle\left(\begin{array}[]{c}q\\ 1\\ \end{array}\right)=\left(\begin{array}[]{cc}R&c\\ 0\cdots 0&1\\ \end{array}\right)\left(\begin{array}[]{c}p\\ 1\\ \end{array}\right).$

(12)

Let

$$A=\left(\begin{array}[]{cc}R&c\\ 0\cdots 0&1\\ \end{array}\right).$$

Trivially, $A$ is linear a linear operator. The matrix $A$ is the most often used form for recording the affine registration.

Let ${\bf p}_{i}$ be the $i$-th landmark and its corresponding affine transformed points ${\bf q}_{i}$. Then we rewrite (12) as

$\displaystyle\underbrace{\left(\begin{array}[]{ccc}{\bf q}_{1}&\cdots&{\bf q}_{n}\end{array}\right)}_{Q}=\left(\begin{array}[]{cc}R&c\\ \end{array}\right)\underbrace{\left(\begin{array}[]{ccc}{\bf p}_{1}&\cdots&{\bf p}_{n}\\ 1&\cdots&1\end{array}\right)}_{P}.$

Then the least squares estimation is given as

$$\left(\begin{array}[]{cc}\widehat{R}&\widehat{c}\\ \end{array}\right)=QP^{\prime}(PP^{\prime})^{-1}.$$

Then the points ${\bf p}_{i}$ are mapped to $\widehat{R}{\bf p}_{i}+\widehat{c}$, which may not coincide with ${\bf q}_{i}$ in general. In practice, landmarks are automatically identified from T1-weighted MRI.

Using SPHARM, it is possible to approximately register surfaces with different mesh topology without any optimization. The crude alignment can be done by coinciding the first order ellipsoid meridian and equator in the SPHARM representation (Gerig et al., 2001; Styner et al., 2006). However, this can be improved. Consider SPHARM representation of surface $H=(h_{1},h_{2},h_{3})$ with spherical angles $\Omega$ given by

$$h_{i}(\Omega)=\sum_{l=0}^{k}\sum_{m=-l}^{l}h_{lm}^{i}Y_{lm}(\Omega),$$

where $(v_{1},v_{2},v_{3})$ are the coordinates of mesh vertices and SPHARM coefficient $h_{lm}^{i}=\langle v_{i},Y_{lm}\rangle$. Consider another SPHARM representation $J=(j_{1},j_{2},j_{3})$ obtained from mesh coordinate $w_{i}$:

$$j_{i}(\Omega)=\sum_{l=0}^{k}\sum_{m=-l}^{l}w_{lm}^{i}Y_{lm}(\Omega),$$

where $w_{lm}^{i}=\langle w_{i},Y_{lm}\rangle$. Suppose the surface $h_{i}$ is deformed to $h_{i}+d_{i}$ by the amount of displacement $d_{i}$. We wish to find $d_{i}$ that minimizes the discrepancy between $h_{i}+d_{i}$ and $j_{i}$ in the subspace $\mathcal{H}_{k}$ spanned using up to the $k$-th degree spherical harmonics. It can be shown that

$\displaystyle\arg\min_{d_{i}\in\mathcal{H}_{k}}\Big{|}\Big{|}\widehat{h_{i}}+d_{i}-\widehat{j_{i}}\Big{|}\Big{|}=\sum_{l=0}^{k}\sum_{m=-l}^{l}(w_{lm}^{i}-v_{lm}^{i})Y_{lm}(\Omega).$

(14)

This implies that the optimal displacement of matching two surfaces is obtained by simply taking the difference between two SPHARM and matching coefficients of the same degree and order. Then a specific point $\widehat{h_{i}}(\Omega)$ in one surface corresponds to $\widehat{j_{i}}(\Omega)$ in the other surface. We refer to this point-to-point surface matching as the SPHARM correspondence (Chung et al., 2008). Unlike other surface registration methods (Chung et al., 2005; Robbins, 2003; Thompson & Toga, 1996), it is not necessary to consider an additional cost function that guarantees the smoothness of the displacement field since the displacement field $d=(d_{1},d_{2},d_{3})$ is already a linear combination of smooth basis functions.

Diffeomorphic image registration is a recently popular technique for registering volume and surface data (Avants et al., 2008; Miller & Qiu, 2009; Vaillant et al., 2007; Qiu & Miller, 2008; Yang et al., 2011). From the affine transformed individual surfaces $\mathcal{M}_{j}$, an additional nonlinear surface registration to the template using the large deformation diffeomorphic metric mapping (LDDMM) framework can be performed (Miller & Qiu, 2009; Vaillant et al., 2007; Qiu & Miller, 2008; Yang et al., 2011). In LDDMM, given the template surface $\mathcal{M}$, the diffeomorphic transformations, which are one-to-one, smooth forward, and inverse transformation, are constructed as follows. We estimate the diffeomorphism $d$ between them as a smooth streamline given by the Lagrangian evolution:

$$\frac{\partial d}{\partial t}(x,t)=v\circ d(x,t)$$

with $d(x,0)=x,\;t\in[0,1]$ for time-dependent velocity field $v$. Note the surfaces $\mathcal{M}_{j}$ and $\mathcal{M}$ are the start and end points of the diffeomorphism, i.e. $\mathcal{M}_{j}\circ d(\cdot,0)=\mathcal{M}_{j}$ and $\mathcal{M}_{j}\circ d(\cdot,1)=\mathcal{M}$. By solving the evolution equation numerically, we obtain the diffeomorphism. By averaging the inverse deformation fields from the template $\mathcal{M}$ to individual subjects, we may obtain yet another final more refined template. The vector fields $v$ are constrained to be sufficiently smooth to generate diffeomorphic transformations over finite time (Dupuis et al., 1998). Figure 2 shows the resulting displacement vector fields of warping from the template to two hippocampal surfaces. In the deformation-based morphometry, variability in the displacement is used to characterize surface growth and differences (Chung, 2012).

Once we extract both the outer and inner cortical surface meshes, cortical thickness can be computed at each mesh vertex. The cortical thickness is defined as the distance between the corresponding vertices between the inner and outer surfaces (MacDonald et al., 2000). There are many different computational approaches in measuring cortical thickness. In one approach, the vertices on the inner triangular mesh are deformed to fit the outer surface by minimizing a cost function that involves bending, stretch and other topological constraints (Chung, 2012). There is also an alternate method for automatically measuring cortical thickness based on the Laplace equation (Jones et al., 2000).

The average cortical thickness for each individual is about 3 mm (Henery & Mayhew, 1989). Cortical thickness varies from 1 to 4 mm depending on the location of the cortex. In normal brain development, it is highly likely that the change of cortical thickness may not be uniform across the cortex. Since different clinical populations are expected to show different patterns of cortical thickness variations, cortical thickness has also been used as a quantitative index for characterizing a clinical population (Chung et al., 2005). Cortical thickness varies locally by region and is likely to be influenced by aging, development and disease (Barta et al., 2005). By analyzing how cortical thickness varies between clinical and non-clinical populations, we can locate regions on the brain related to a specific pathology.

As in the case of local volume changes in the deformation-based morphometry, the rate of cortical surface area expansion or reduction may not be uniform across the cortical surface (Chung, 2012). Suppose that cortical surface $\mathcal{M}$ is parameterized by the parameters ${\bf u}=(u^{1},u^{2})$ such that any point ${\bf x}\in\mathcal{M}$ can be written as ${\bf x}={\bf X}({\bf u})$. Let ${\bf X}_{i}=\partial{\bf X}/\partial u^{i}$ be the partial derivative vectors. The Riemannian metric tensor $g_{ij}$ is then given by the inner product between two vectors ${\bf X}_{i}$ and ${\bf X}_{j}$, i.e.

$$g_{ij}(t)=\langle{\bf X}_{i},{\bf X}_{j}\rangle.$$

The tensor $g_{ij}$ measures the amount of the deviation of the cortical surface from a flat Euclidean plane and can be used to measure lengths, angles and areas on the cortical surface.

Let $g=(g_{ij})$ be a $2\times 2$ matrix of metric tensors. The total surface area of the cortex $\mathcal{M}$ is then given by

$$\int_{D}\sqrt{\det g}\;d{\bf u},$$

where $D=X^{-1}(\mathcal{M})$ is the parameter space (Kreyszig, 1959). The term $\sqrt{\det g}$ is called the surface area element and it measures the transformed area of the unit square in the parameter space $D$ via transformation $X:D\to\mathcal{M}$. The surface area element can be considered as the generalization of the Jacobian determinant, which is used in measuring local volume changes in the tensor-based morphometry (Chung, 2012).

Instead of using the metric tensor formulation, it is possible to quantify local surface area change in terms of the areas of the corresponding triangles in surface meshes. However, this formulation assigns the computed surface area to each face instead of each vertex. This causes a problem in both surface-based smoothing and statistical analysis, where values are required to be given on vertices. Interpolating scalar values on vertices from face values can be done by the weighted average of face values. It is not hard to develop surface-based smoothing and statistical analysis on face values, as a form of dual formulation, but the cortical thickness and the curvature metric are defined on vertices so we end up with two separate approaches: one for metrics defined on vertices and the other for metrics defined on faces. Therefore, the metric tensor approach provides a better unified quantitative framework for the subsequent statistical analysis.

The principal curvatures characterize the shape and location of the sulci and gyri, which are the valleys and crests of the cortical surfaces (Bartesaghi & Sapiro, 2001; Joshi et al., 1995; Khaneja et al., 1998). By measuring the curvature changes, rapidly folding and cortical regions can be localized. The principal curvatures $\kappa_{1}$ and $\kappa_{2}$ can be represented as functions of $\beta_{i}$ in quadratic surface (1) (Boothby, 1986; Kreyszig, 1959).

Local volume can be computed using the determinant of the Jacobian of deformation and used in detecting the regions of brain tissue growth and loss in brain development (Chung, 2012). Compared to the local surface area change, the local volume change measurement is more sensitive to small deformation of the brain. If a unit cube increases its sides by one, the surface area will increase by $2^{2}-1=3$ while the volume will increase by $2^{3}-1=7$. Therefore, the statistical analysis based on the local volume change is at least twice more sensitive compared to that of the local surface area change. So the local volume change should be able to pick out gray matter tissue growth patterns even when the local surface area change may not.

The gray matter can be considered as a thin shell bounded by two surfaces with varying cortical thickness. In most deformable surface algorithms like FreeSurfer, each triangle on the outer surface has a corresponding triangle on the inner surface. Let ${\bf p}_{1},{\bf p}_{2},{\bf p}_{3}$ be the three vertices of a triangle on the outer surface and ${\bf q}_{1},{\bf q}_{2},{\bf q}_{3}$ be the corresponding three vertices on the inner surface such that ${\bf p}_{i}$ is linked to ${\bf q}_{i}$. Then the volume of the triangular prism is given by the sum of the determinants

$$D({\bf p}_{1},{\bf p}_{2},{\bf p}_{3},{\bf q}_{1})+D({\bf p}_{2},{\bf p}_{3},{\bf q}_{1},{\bf q}_{2})+D({\bf p}_{3},{\bf q}_{1},{\bf q}_{2},{\bf q}_{3})$$

where

$$D({\bf a},{\bf b},{\bf c},{\bf d})=|\det({\bf a}-{\bf d},{\bf b}-{\bf d},{\bf c}-{\bf d})|/6$$

is the volume of a tetrahedron whose vertices are $\{{\bf a},{\bf b},{\bf c},{\bf d}\}$. Afterwards, the total gray matter volume can be estimated by summing the volumes of all triangular prisms (Chung, 2012).

By formulating Gaussian kernel smoothing as a solution of a diffusion equation on a Riemannian manifold, the Gaussian kernel smoothing approach can be generalized to an arbitrary curved surface. This generalization is called diffusion smoothing and was first introduced in the analysis of fMRI data on the cortical surface (Andrade et al., 2001) and cortical thickness (Chung et al., 2001) in 2001.

It can be shown that Gaussian kernel smoothing (15) is the integral solution of the $n$-dimensional diffusion equation

$\displaystyle\frac{\partial F}{\partial t}=\Delta F$

(16)

with the initial condition $F({\bf x},0)=f({\bf x})$, where

$$\Delta=\frac{\partial^{2}}{\partial x_{1}^{2}}+\cdots+\frac{\partial^{2}}{\partial x_{n}^{2}}$$

is the Laplacian in $n$-dimensional Euclidean space. Hence the Gaussian kernel smoothing is equivalent to the diffusion of the initial data $f({\bf x})$ after time $t$.

Diffusion equations have been widely used in image processing as a form of noise reduction starting with Perona and Malik in 1990 (Perona & Malik, 1990). Numerous diffusion techniques have been developed for surface data smoothing (Andrade et al., 2001; Chung, 2012; Cachia, Mangin, Riviére, Kherif, Boddaert, Andrade, Papadopoulos-Orfanos, Poline, Bloch, Zilbovicius, Sonigo, Brunelle & Régis, 2003; Cachia, Mangin, Riviére, Papadopoulos-Orfanos, Kherif, Bloch & Régis, 2003; Chung et al., 2005; Joshi et al., 2009). When applying diffusion smoothing on curved surfaces, the smoothing somehow has to incorporate the geometrical features of the curved surface and the Laplacian $\Delta$ should change accordingly. The extension of the Euclidean Laplacian to an arbitrary Riemannian manifold is called the Laplace-Beltrami operator (Arfken, 2000; Kreyszig, 1959). The approach taken in (Andrade et al., 2001) is based on a local flattening of the cortical surface and estimating the planar Laplacian, which may not be as accurate as the cotan estimation based on the finite element method (FEM) given in (Chung et al., 2001). Further, the FEM approach completely avoids the use of surface flattening and parameterization; thus, it is more robust.

For given Riemannian metric tensor $g_{ij}$, the Laplace-Beltrami operator $\Delta$ is defined as

$\displaystyle\Delta F=\sum_{i,j}\frac{1}{|g|^{1/2}}\frac{\partial}{\partial u^{i}}\Big{(}|g|^{1/2}g^{ij}\frac{\partial F}{\partial u^{j}}\Big{)},$

(17)

where $(g^{ij})=g^{-1}$ (Arfken, 2000). Note that when $g$ becomes an identity matrix, the Laplace-Beltrami operator reduces to the standard 2D Laplacian:

$$\Delta F=\frac{\partial^{2}F}{\partial(u^{1})^{2}}+\frac{\partial^{2}F}{\partial(u^{2})^{2}}.$$

Using FEM on the triangular cortical mesh, it is possible to estimate the Laplace-Beltrami operator as the linear weights of neighboring vertices using the cotan formulation, which is first given in (Chung et al., 2001).

Let ${\bf p}_{1},\cdots,{\bf p}_{m}$ be $m$ neighboring vertices around the central vertex ${\bf p}={\bf p}_{0}$ (Figure 3). Then the estimated Laplace-Beltrami operator is given by

$\displaystyle\widehat{\Delta F}({\bf p})=\sum_{i=1}^{m}w_{i}\big{(}F({\bf p}_{i})-F({\bf p})\big{)}$

(18)

with the weights

$$w_{i}=\frac{\cot\theta_{i}+\cot\phi_{i}}{\sum_{i=1}^{m}\|T_{i}\|},$$

where $\theta_{i}$ and $\phi_{i}$ are the two angles opposite to the edge connecting ${\bf p}_{i}$ and ${\bf p}$, and $\|T_{i}\|$ is the area of the $i$-th triangle (Figure 3).

FEM estimation (18) is an improved formulation from the previous attempt in diffusion smoothing (Andrade et al., 2001), where the Laplacian is simply estimated as the planar Laplacian after locally flattening the triangular mesh consisting of nodes ${\bf p}_{0},\cdots,{\bf p}_{m}$ onto a flat plane. Afterwards, the finite difference (FD) scheme can be used to iteratively solve the diffusion equation at each vertex ${\bf p}$:

$$\frac{F({\bf p},t_{n+1})-F({\bf p},t_{n})}{t_{n+1}-t_{n}}=\widehat{\Delta}F({\bf p},t_{n}),$$

with the initial condition $F({\bf p},0)=f({\bf p})$ and fixed $\delta t=t_{n+1}-t_{n}$. After $N$ iterations, the FD gives the diffusion of the initial data $f$ after time $N\delta t$. If the diffusion were applied to Euclidean space, it would be approximately equivalent to Gaussian kernel smoothing with

$$\mbox{FWHM}=4(\ln 2)^{1/2}\sqrt{N\delta t}.$$

For large meshes, computing the linear weights for the Laplace-Beltrami operator takes a fair amount of time, but once the weights are computed, it is applied repeatedly throughout the iterations as a matrix multiplication. Unlike Gaussian kernel smoothing, diffusion smoothing is an iterative procedure.

Diffusion smoothing use FEM and FD, which are known to suffer numerical instability if sufficiently small step size is not chosen in the forward Euler scheme. To remedy the problem associated with diffusion smoothing, iterative kernel smoothing³³3MATLAB package: www.stat.wisc.edu/$\sim$mchung/softwares/hk/hk.html. was introduced (Chung et al., 2005). The method has been used in smoothing various cortical surface data: cortical curvatures (Luders, Thompson, Narr, Toga, Jancke & Gaser, 2006; Gaser et al., 2006), cortical thickness (Luders, Narr, Thompson, Rex, Woods, DeLuca, Jancke & Toga, 2006; Bernal-Rusiel et al., 2008), hippocampus (Shen et al., 2006; Zhu et al., 2007), magnetoencephalography (MEG) (Han et al., 2007) and functional-MRI (Hagler Jr. et al., 2006; Jo et al., 2007). This and its variations are probably the most widely used method for smoothing brain surface data at this moment. In iterated kernel smoothing, kernel weights are spatially adapted to follow the shape of the heat kernel in a discrete fashion.

The $n$-th iterated kernel smoothing of signal $f\in L^{2}(\mathcal{M})$ with kernel $K_{\sigma}$ is defined as

$$K_{\sigma}^{(n)}*f(p)=\underbrace{K_{\sigma}*\cdots*K_{\sigma}}_{n\mbox{ \small{times} }}*f(p),$$

where $\sigma$ is the bandwidth of the kernel. If $K_{\sigma}$ is a heat kernel, we have the following iterative relation (Chung et al., 2005):

$\displaystyle K_{\sigma}*f(p)=K_{\sigma/n}^{(n)}*f(p).$

(19)

The relation (19) shows that kernel smoothing with large bandwidth $\sigma$ can be decomposed into $n$ repeated applications of kernel smoothing with smaller bandwidth $\sigma/n$. This idea can be used to approximate the heat kernel. When the bandwidth is small, the heat kernel behaves like the Dirac-delta function and, using the parametrix expansion (Rosenberg, 1997; Wang, 1997), we can approximate it locally using the Gaussian kernel.

Let ${\bf p}_{1},\cdots,{\bf p}_{m}$ be $m$ neighboring vertices of vertex ${\bf p}={\bf p}_{0}$ in the mesh. The geodesic distance $d({\bf p},{\bf p}^{\prime})$ between ${\bf p}$ and its adjacent vertex ${\bf p}_{i}$ is the length of edge between these two vertices in the mesh. Then the discretized and normalized heat kernel is given by

$$W_{\sigma}({\bf p},{\bf p}_{i})=\frac{\exp\big{(}-\frac{d({\bf p},{\bf p}_{i})^{2}}{4\sigma}\big{)}}{\sum_{j=0}^{m}\exp\big{(}-\frac{d({\bf p}-{\bf p}_{j})^{2}}{4\sigma}\big{)}}.$$

Note $\sum_{i=0}^{m}W_{\sigma}({\bf p},{\bf p}_{i})=1$. For small bandwidth, all the kernel weights are concentrated near the center, so we only need to worry about the first neighbors of a given vertex in a surface mesh. The discrete version of heat kernel smoothing on a triangle mesh is then given by

$$W_{\sigma}*f(p)=\sum_{i=0}^{m}W_{\sigma}(p,p_{i})f(p_{i}).$$

The discrete kernel smoothing should converge to heat kernel smoothing as the mesh resolution increases. This is the form of the Nadaraya-Watson estimator (Chaudhuri & Marron, 2000) in statistics. Instead of performing a single kernel smoothing with large bandwidth $n\sigma$, we perform $n$ iterated kernel smoothing with small bandwidth $\sigma$ as follows $W_{\sigma}^{(n)}*f(p)$.

The recently proposed heat kernel smoothing⁴⁴4MATLAB package: brainimaging.waisman.wisc.edu/$\sim$chung/lb. framework constructs the heat kernel analytically using the eigenfunctions of the Laplace-Beltrami operator (Seo et al., 2010). This method avoids the need for the linear approximation used in iterative kernel smoothing that compounds the approximation error at each iteration. The method represents isotropic heat diffusion analytically as a series expansion so it avoids the numerical convergence issues associated with solving the diffusion equations numerically (Andrade et al., 2001; Chung, 2012; Joshi et al., 2009). This technique is different from other diffusion-based smoothing methods in that it bypasses the various numerical problems such as numerical instability, slow convergence, and accumulated linearization error.

Although recently there have been a few studies that introduce heat kernel in computer vision and machine learning (Belkin et al., 2006), they mainly use heat kernel to compute shape descriptors (Sun et al., 2009; Bronstein & Kokkinos, 2010) or to define a multi-scale metric (de Goes et al., 2008). These studies did not use heat kernel in smoothing functional data on manifolds. Further, most kernel methods in machine learning deal with the linear combination of kernels as a solution to penalized regressions, which significantly differs from the heat kernel smoothing framework, which does not have a penalized cost function. There are log-Euclidean and exponential map frameworks on manifolds, where the main interest is in computing the Fréchet mean along the tangent space (Davis et al., 2010; Fletcher et al., 2004). Such approaches or related methods in (Gerber et al., 2009), the Nadaya-Watson type of kernel regression is reformulated to learn the shape or image means in a population. In the heat kernel smoothing framework, we are not dealing with manifold data but scalar data defined on a manifold, so there is no need for exploiting the manifold structure of the data itself.

Let $\Delta$ be the Laplace-Beltrami operator on $\mathcal{M}$. Solving the eigenvalue equation

$\displaystyle\Delta\psi_{j}=-\lambda\psi_{j},$

(20)

we order eigenvalues

$$0=\lambda_{0}<\lambda_{1}\leq\lambda_{2}\leq\cdots,$$

and corresponding eigenfunctions $\psi_{0},\psi_{1},\psi_{2},\cdots$ (Rosenberg, 1997; Chung et al., 2005; Lévy, 2006; Shi et al., 2009). The eigenfunctions $\psi_{j}$ form an orthonormal basis in $L^{2}(\mathcal{M})$, the space of square integrable functions in $\mathcal{M}$. Figure 4 shows the first four LB-eigenfunctions on a hippocampal surface.

There is extensive literature on the use of eigenvalues and eigenfunctions of the Laplace-Beltrami operator in medical imaging and computer vision (Lévy, 2006; Qiu et al., 2006; Reuter et al., 2009; Reuter, 2010; Zhang et al., 2007, 2010). The eigenvalues have been used in caudate shape discriminators (Niethammer et al., 2007). Qiu et al. used eigenfunctions in constructing splines on cortical surfaces (Qiu et al., 2006). Reuter used the topological features of eigenfunctions (Reuter, 2010). Shi et al. used the Reeb graph of the second eigenfunction in shape characterization and landmark detection in cortical and subcortical structures (Shi et al., 2008, 2009). Lai et al. used the critical points of the second eigenfunction as anatomical landmarks for colon surfaces (Lai et al., 2010).

Using the eigenfunctions, heat kernel $K_{\sigma}(p,q)$ is defined as

$$K_{\sigma}(p,q)=\sum_{j=0}^{\infty}e^{-\lambda_{j}\sigma}\psi_{j}(p)\psi_{j}(q),$$

(21)

where $\sigma$ is the bandwidth of the kernel. The heat kernel is the generalization of a Gaussian kernel. Then heat kernel smoothing of functional measurement $Y$ is defined as

$$K_{\sigma}\ast Y(p)=\sum_{j=0}^{\infty}e^{-\lambda_{j}\sigma}\beta_{j}\psi_{j}(p),$$

(22)

where $\beta_{j}=\langle Y,\psi_{j}\rangle$ are Fourier coefficients (Chung et al., 2005). Kernel smoothing $K_{\sigma}*Y$ is taken as the estimate for the unknown mean signal $\theta$. The degree for truncating the series expansion can be automatically determined using the forward model selection procedure (Chung et al., 2008). Figure 4 shows the heat kernel smoothing results with the bandwidth $\sigma=0.5$ and $k=500$ number of LB-eigenfunctions.

Unlike previous approaches to heat diffusion (Andrade et al., 2001; Chung, 2012; Joshi et al., 2009; Tasdizen et al., 2006), heat kernel smoothing avoids the direct numerical discretization of the diffusion equation. Instead, we discretize the basis functions of given manifold $\mathcal{M}$ by solving for the eigensystem (20) and obtain $\lambda_{j}$ and $\psi_{j}$. This provides more robust stable smoothing results compared to diffusion smoothing or iterated kernel smoothing approaches.

We set up a GLM at each mesh vertex. Let $y_{i}$ be the response variable, which is mainly coming from images and ${\bf x}_{i}=(x_{i1},\cdots,x_{ip})$ to be the variables of interest and ${\bf z}_{i}=(z_{i1},\cdots,z_{ik})$ to be nuisance variables corresponding to the $i$-th subject. Assume there are $n$ subjects, i.e., $i=1,\cdots,n$. We are interested in testing the significance of variables ${\bf x}_{i}$ while accounting for nuisance covariates ${\bf z}_{i}$. Then we set up GLM

$$y_{i}={\bf z}_{i}{\boldsymbol{\lambda}}+{\bf x}_{i}{\boldsymbol{\beta}}+\epsilon_{i}$$

where ${\boldsymbol{\lambda}}=(\lambda_{1},\cdots,\lambda_{k})^{\prime}$ and ${\boldsymbol{\beta}}=(\beta_{1},\cdots,\beta_{p})^{\prime}$ are unknown parameter vectors to be estimated. We assume $\epsilon$ to be the usual zero mean Gaussian noise, although the distributional assumption is not required for the least squares estimation. We test hypotheses

$$H_{0}:{\boldsymbol{\beta}}=0\mbox{ vs. }H_{1}:{\boldsymbol{\beta}}\neq 0.$$

Subsequently the inference is done by constructing the $F$-statistic with $p$ and $n-p-k$ degrees of freedom. GLMs have been used in quantifying cortical thickness, for instance, in child development (Chung, 2012; Chung et al., 2005) and amygdala shape differences in autism (Chung, 2012).

In the hippocampus case study, the first T1-weighted MRI scans are taken at $11.6\pm 3.7$ years for $n=124$ children using a 3T GE SIGNA scanner. Variables age and gender are available. We also have variable income, which is a binary dummy variable indicating whether the subjects are from high- or low-income families. A total of 124 children and adolescents are from high- ($>75000\$$; $n=86$) and low-income ($<35000\$$, $n=38$) parents respectively. In addition to this cross-sectional data, longitudinal data was available for 82 of these subjects ($n=66$, $>75000\$$; $n=16$, $<35000\$$). The second MRI scan was acquired for these 82 subjects about 2 years later at $14\pm 3.9$ years. For now, we will simply ignore the correlation between the scans within a subject and will treat them as independent.

On the template surface, we have the displacement vector fields of mapping from the template to individual subjects (Figure 2). We take the length of the surface displacement, denoted as deformation, with respect to the template as the response variable. The displacement measures the shape difference with respect to the template. However, since the length measurement is noisy, surface-based smoothing is necessary. We have used heat kernel smoothing to smooth out noise with the bandwidth 1 and 500 LB-eigenfunctions. Then we set up the GLM:

$${\tt deformation}=\lambda_{1}+\lambda_{2}{\tt age}+\beta_{1}{\tt income}+\epsilon$$

and test for the significance of $\beta_{1}$ at each mesh vertex. Figure 5-left shows the $F$-statistic result on testing $\beta_{1}$.

The multivariate general linear models (MGLMs) have been also used in modeling multivariate imaging features on brain surfaces. These models generalize a widely used univariate GLM by incorporating vector valued responses and explanatory variables (Anderson, 1984; Friston et al., 1995; Worsley et al., 1996, 2004; Taylor & Worsley, 2008; Chung, 2012). Hotelling’s $T^{2}$ statistic is a special case of MGLM that has been used primarily for inference on surface shapes, deformations and multi-scale features (Cao & Worsley, 1999; Chung, 2012; Gaser et al., 1999; Joshi, 1998; Kim et al., 2012; Thompson et al., 1997). An example of this approach is Hotelling’s $T^{2}$-statistic applied in determining the 3D brain morphology of HIV/AIDS patients (Lepore et al., 2006). Hotelling’s $T^{2}$-statistic is also applied to 2D deformation tensor at each mesh vertex on the hippocampal surface as a way to characterize AlzheimerÕs diseases (Wang et al., 2011).

Suppose there are a total of $n$ subjects and $p$ multivariate features of interest at each voxel. For MGLM to work, $n$ should be significantly larger than $p$. Let ${\bf J}_{n\times p}=(J_{ij})$ be the measurement matrix, where $J_{ij}$ is the measurement for subject $i$ and for the $j$-th feature. The subscripts denote the dimension of the matrix. All the measurements over subjects for the $j$-th feature are denoted as ${\bf x}_{j}=(J_{1j},\cdots,J_{nj})^{\prime}$. The measurement vector for the $i$-th subject is denoted as ${\bf y}_{i}=(J_{i1},\cdots,J_{ip})^{\prime}$. ${\bf y}_{i}$ is expected to be distributed identically and independently over subjects. Note that

$${\bf J}=({\bf x}_{1},\cdots,{\bf x}_{p})=({\bf y}_{1},\cdots,{\bf y}_{n})^{\prime}.$$

We may assume the covariance matrix of ${\bf y}_{i}$ to be

$$\mbox{Cov}({\bf y}_{1})=\cdots=\mbox{Cov}({\bf y}_{n})={\bf\Sigma}_{p\times p}=(\sigma_{kl}).$$

With these notations, we now set up the following MGLM at each mesh vertex:

$\displaystyle{\bf J}_{n\times p}={\bf X}_{n\times k}{\bf B}_{k\times p}+{\bf Z}_{n\times q}{\bf G}_{q\times p}+{\bf U}_{n\times p}{\bf\Sigma}_{p\times p}^{1/2}.$

(24)

$\bf X$ is the matrix of contrasted explanatory variables while $\bf B$ is the matrix of unknown coefficients. Nuisance covariates are in matrix $\bf Z$ and the corresponding coefficients are in matrix $\bf G$. The components of Gaussian random matrix $\bf U$ are independently distributed with zero mean and unit variance. Symmetric matrix $\bf\Sigma^{1/2}$ is the square root of the covariance matrix accounting for the spatial dependency across different voxels. In MGLM (24), we are usually interested in testing hypotheses

$$H_{0}:{\bf B}=0.\;\mbox{\em vs.}\;H_{1}:{\bf B}\neq 0.$$

The parameter matrices in the model are estimated via the least squares method. The multivariate test statistics such as Lawley-Hotelling trace or Roy’s maximum root are used to test the significance of ${\bf B}$. When there is only one voxel, i.e. $p=1$, these multivariate test statistics collapse to Hotelling’s $T^{2}$ statistic (Worsley et al., 2004; Chung, 2012).

GLM are usually fitted in each voxel separately. Instead of fitting GLM at each voxel, one may be tempted to fit the model in the whole brain surface. For FreeSurfer meshes, we need to fit GLM over 300000 vertices, which causes the small-$n$ large-$p$ problem (Friston et al., 1995; Schäfer & Strimmer, 2005; Chung, 2013).

Let ${\bf y}_{j}$ be the measurement vector at the $j$-th vertex. Assume there are $n$ subjects and total $p$ vertices in the surface. We have the same design matrix ${\bf Z}$ for all $p$ vertices. Then we need to estimate the parameter vector ${\boldsymbol{\lambda}}_{j}$ in

$\displaystyle{\bf y}_{j}={\bf Z}{\boldsymbol{\lambda}_{j}}$

(25)

for each $j$. Instead of solving (25) separately at each vertex, we combine all of them together so that we have matrix equation

$\displaystyle\underbrace{[{\bf y}_{1},\cdots,{\bf y}_{m}]}_{\bf Y}={\bf Z}\underbrace{[{\boldsymbol{\lambda}}_{1},\cdots,{\boldsymbol{\lambda}}_{m}]}_{\boldsymbol{\Lambda}}.$

(26)

The least squares estimation of the parameter matrix ${\boldsymbol{\Lambda}}$ is given by

$$\widehat{\boldsymbol{\Lambda}}=({\bf Z}^{\prime}{\bf Z})^{-1}{\bf Z}^{\prime}{\bf Y}.$$

Note that ${\bf Z}$ is of size $n$ by $p$ and ${\bf Z}^{\prime}{\bf Z}$ is only invertible when $n\ll p$. The least squares estimation does not provide robust parameter estimates for $n\ll p$, which is the usual case in surface modeling. For small-$n$ large-$p$ problem, GLM need to be regularized using the L1-norm penalty (Banerjee et al., 2006, 2008; Chung, 2013; Friedman et al., 2008; Huang et al., 2010; Mazumder & Hastie, 2012).

So far we have only dealt with an imaging data set where the parameters of the model are fixed and do not vary across subjects and scans. Such fixed-effects models are inadequate in modeling the within-subject dependency in longitudinally collected imaging data. However, mixed-effects models can explicitly model such dependency (Fox, 2002; Milliken & Edland, 2000; Molenberghs & Verbeke, 2005; Pinehiro & Bates, 2002). There are three advantages of the mixed-effects model over the usual fixed-effects model. It explicitly models individual growth patterns and accommodates an unequal number of follow-up image scans per subject and unequal time intervals between scans.

The longitudinal outcome $Y_{i}$ from the $i$-th subject is modeled using the mixed-effects model (Milliken & Edland, 2000) as

$\displaystyle Y_{i}=X_{i}\beta+Z_{i}\gamma_{i}+e_{i},$

(27)

where $\beta$ is the fixed-effects shared by all subjects. $\gamma_{i}$ is the subject-specific random-effects and $e_{i}\sim N(0,\sigma^{2})$ is independent and identically distributed noise. $X_{i}$ and $Z_{i}$ are the design matrices corresponding to the fixed and random effects respectively for the $i$-th subject. We assume $\gamma_{i}\sim N(0,\Gamma)$ and $\epsilon_{i}\sim N(0,\Sigma_{i})$ with some covariance matrices $\Gamma$ and $\Sigma_{i}$. Hierarchically we are modeling (27) as

$$Y_{i}|\gamma_{i}\sim N(X_{i}\beta+Z_{i}\gamma_{i},\Sigma_{i}),\;\gamma_{i}\sim N(0,\Gamma).$$

$\Gamma$ accounts for covariance among random effect terms. The within-subject variability between the scans is expected to be smaller than between-subject variability and explicitly modeled by $\Sigma_{i}$. The covariance of $\gamma_{i}$ and $\epsilon$ are expected to have block diagonal structure such that there is no correlation among the scans of different subjects while there is high correlation between the scans of the same subject:

$$\mathbb{V}\left(\begin{array}[]{c}\gamma_{i}\\ \epsilon_{i}\\ \end{array}\right)=\left(\begin{array}[]{ccc}\Gamma&0\\ 0&\Sigma_{i}\\ \end{array}\right).$$

Subsequently, the overall covariance of $Y_{i}$ is given by

$$\mathbb{V}Y_{i}=Z_{i}\Gamma Z_{i}^{\prime}+\Sigma_{i}.$$

The random-effect contribution is $Z_{i}\Gamma Z_{i}^{\prime}$ while the within-subject contribution is $\Sigma_{i}$.

The parameters and the covariance matrices can be estimated via the restricted maximum likelihood (REML) method (Fox, 2002; Pinehiro & Bates, 2002). The most widely used tools for fitting the mixed-effects model are the nlme library in the R statistical package (Pinehiro & Bates, 2002). However, there is no need to use R to fit the mixed-effects model. Keith Worsley has implemented the REML procedure in the SurfStat package⁶⁶6The MATLAB package is available at www.math.mcgill.ca/keith/surfstat. (Chung, 2012; Worsley et al., 2009).

Here we briefly explain how to set up a longitudinal mixed-effect model in practice. In the usual fixed-effect model, we have a linear model containing the fixed-effect term ${\tt age}_{i}$ for the $i$-th subject:

$\displaystyle y_{i}=\beta_{0}+\beta_{1}{\tt age}_{i}+\epsilon_{i},$

(28)

where $\epsilon_{i}$ is assumed to follow independent Gaussian. In (28), every subject has identical growth trajectory $\beta_{0}+\beta_{1}{\tt age}$, which is unrealistic. Biologically, each subject is expected to have its own unique growth trajectory. So we assume each subject to have its own intercept $\beta_{0}+\gamma_{i0}$ and slope $\beta_{1}+\gamma_{i1}$:

$\displaystyle y_{i}=\beta_{0}+\gamma_{i0}+(\beta_{1}+\gamma_{i1}){\tt age}_{i}+\epsilon_{i}.$

(29)

It is reasonable to assume the random vector $\gamma=(\gamma_{i0},\gamma_{i1})$ to be multivariate normal. The model (29) can be decomposed into fixed- and random-effect terms:

$\displaystyle y_{i}=(\beta_{0}+\beta_{1}{\tt age}_{i})+(\gamma_{i0}+\gamma_{i1}{\tt age}_{i})+\epsilon_{i}.$

(30)

The fixed-effect term $\beta_{0}+\beta_{1}{\tt age}_{i}$ models the linear growth of the population while the random-effect term $\gamma_{i0}+\gamma_{i1}{\tt age}_{i}$ models the subject specific growth variations. Incorporating additional factors and interaction terms are done similarly.

In the hippocampus case study, the first MRI scans are taken at $11.6\pm 3.7$ years while the second scans are taken at $14\pm 3.9$ years. We are interested in determining the effects of income level on the shape of the hippocampus. In Section 6.1, we treated the second scans as if they came from independent subjects and modeled them using the fixed-effects model. Now we explicitly incorporate the dependency of repeated scans of the same subject. It is necessary to explicitly model the within-subject variability that is expected to be smaller than between-subject variability. This can be done by introducing a random-effect term. The resulting F-statistic maps are given in Figure 5-right. However, we did not detect any region that is affected by income. Thus, we tested the age and income interaction and found the regions of strong interaction (Figure 6).

Since we need to set up a GLM on every mesh vertex, it becomes a multiple comparisons problem. Correcting for multiple comparisons is crucial in determining overall statistical significance in correlated test statistics over the whole surface. For surface data, various methods are proposed: Bonferroni correction, random field theory (Worsley, 1994; Worsley et al., 1996), false discovery rates (FDR) (Benjamini & Hochberg, 1995; Benjamini & Yekutieli, 2001; Genovese et al., 2002), and permutation tests (Nichols & Holmes, 2002). Among many techniques, the random field theory is probably the most natural in relation to surface data smoothing since it is able to explicitly control the amount of smoothing.

The generalization of a continuous stochastic process in $\mathbb{R}^{n}$ to a higher dimensional abstract space is called a random field (Adler & Taylor, 2007; Chung, 2012; Dougherty, 1999; Joshi, 1998; Yaglom, 1987). In the random field theory (Worsley, 1994; Worsley et al., 1996), measurement $Y$ at position $x\in\mathcal{M}$ is modeled as

$$Y(x)=\mu(x)+\epsilon(x),$$

where $\mu$ is the unknown signal to be estimated and $\epsilon$ is the measurement error. The measurement error at each fixed $x$ can be modeled as a random variable. Then the collection of random variables $\{\epsilon(x):x\in\mathcal{M}\}$ is called a random field or stochastic process. A measure-theoretic definition is given in (Adler & Taylor, 2007).

Detecting the regions of statistically significant $\lambda(x)$ for some $x\in\mathcal{M}$ can be done via thresholding the maximum of a random field defined on the cortical surface (Worsley et al., 1996, 1999). For instance, $T$ random field on the surface $\mathcal{M}$ is defined as

$$T(x)=\sqrt{n}\frac{M(x)}{S(x)},$$

where $M$ and $S$ are the sample mean and standard deviation over the $n$ subjects. Under the null hypothesis

$$H_{0}:\mu(x)=0\text{ for all }x\in\mathcal{M},$$

$T(x)$ is distributed as a student’s $T$ with $n-1$ degrees of freedom at each voxel $x$. The pvalue of the local maxima of the $T$ field will give a conservative threshold compared to FDR (Worsley et al., 1996).

For sufficiently high threshold $y$, we can show that

$\displaystyle P\Big{(}\max_{x\in\mathcal{M}}T({\bf x})\geq y\Big{)}\approx\sum_{i=0}^{3}\phi_{i}(\mathcal{M})\rho_{i}(y),$

(31)

where $\rho_{i}$ is the $i$-dimensional EC-density and the Minkowski functional $\phi_{i}$ are

$$\phi_{0}(\mathcal{M})=2,\;\phi_{1}(\mathcal{M})=0,\;\phi_{2}(\mathcal{M}=\|\mathcal{M}\|,\;\phi_{3}(\mathcal{M})=0$$

and $\|\mathcal{M}\|$ is the total surface area of $\mathcal{M}$ (Worsley, 1996a). When diffusion or heat kernel smoothing with given FWHM is applied on surface $\mathcal{M}$, the 0-dimensional and 2-dimensional EC-density becomes

$$\rho_{0}(y)=\int_{y}^{\infty}\frac{\Gamma(\frac{n}{2})}{((n-1)\pi)^{1/2}\Gamma(\frac{n-1}{2})}\Big{(}1+\frac{y^{2}}{n-1}\Big{)}^{-n/2}\;dy,$$

$$\rho_{2}(y)=\frac{1}{\text{FWHM}^{2}}\frac{4\ln 2}{(2\pi)^{3/2}}\frac{\Gamma(\frac{n}{2})}{(\frac{n-1}{2})^{1/2}\Gamma(\frac{n-1}{2})}y\Big{(}1+\frac{y^{2}}{n-1}\Big{)}^{-(n-2)/2}.$$

The excursion probability, which is the probability of obtaining false positives for the one-sided alternate hypothesis, is approximated by the following formula:

$$P\Big{(}\max_{x\in\mathcal{M}}T(x)\geq y\Big{)}\approx 2\rho_{0}(y)+\|\mathcal{M}\|\rho_{2}(y).$$

For smoothing cortical thickness (Chung, 2012; Chung et al., 2005), an FWHM of between 20 to 30 mm is recommended. This FWHM reflects the spatial frequency associated with the sulcal pattern. For measurements on smaller subcortical structures such as hippocampus and amygdala, significantly smaller FWHM is recommended. For amygdala and hippocampus, 0.5–1 mm would be sufficient.

In the hippocampus case study, we did not detect any statistically significant group difference at 0.01 level after correcting for multiple comparisons in both the left and right hippocampi. However, we obtained highly focalized regions of group difference in the growth rate, the interaction between income level and age, in the right hippocampus (corrected pvalue = 0.03). The posterior region is enlarging while the midbody and the anterior parts are shrinking in children from low-income families (Figure 6 and 7). On the other hand, the pattern is opposite for children from high-income families. Note that the right hippocampus is involved in the active maintenance of associations with spatial information (Piekema et al., 2006). Future studies investigating the relationship between family socioeconomic status and spatial information processing measures are warranted.