1 Introduction

Average cortical thickness can vary between 2 to 5 mm across a population of healthy subjects as well as across different brain regions in an individual [1]. Thickness of the human cerebral cortex is an important phenotypical feature and a biomarker for a range of neurological diseases and conditions, and brain development. To perform cortical thickness studies in a large population an automated approach to thickness measurement from T1-weighted MRI scans is essential. Several approaches for computation of cortical thickness are based on first estimating inner gray/white and pial surfaces and then defining the cortical thicknesses based on the distance between the two [2, 3]. The Linked Distance method (LD) in BrainSuite uses distance between corresponding nodes in the two surfaces as a thickness measure. FreeSurfer’s cortical thickness [4] is defined as the average of the shortest distance between the two surfaces computed in both directions. Cortical Pattern Matching (CPM) [5] finds the shortest distance between the inner and pial surfaces using the Eikonal equation. On the other hand, voxel-based methods compute thickness based on line integrals [6, 7], the Laplace equation [1, 8, 9], or using image registration [10]. The accuracy of these approaches is impacted when individual voxels are composed of a mixture of multiple tissue types leading to partial volume effects. The convoluted geometry of the cortex together with the point spread function associated with finite resolution makes partial volume effects inevitable. Despite this, most methods use crisp definitions of cortical boundaries for surface-based calculation. In some cases partial volume effects have been accounted for by modifying the cortical surface boundaries using Eulerian or Lagrangian PDEs [9], registration of inner and pial cortical surfaces [10], closest point distances between inner and pial cortical surfaces measured in both directions [11] and using an electric field model together with a topology preserving level set approach [12]. While these methods account for partial volume effects in defining inner and pial cortical surfaces, they do not explicit account for the actual partial volume fractions that lie between the two boundaries once they are defined.

Here we describe a new thickness calculation method that explicitly models partial volume effects by modifying the Laplace equation (LE) based method of Jones et al. [8]. Rather than using the isotropic LE, we instead use the anisotropic version in which the diffusion coefficient is varied spatially in proportion to the fraction of gray matter in each voxel. Further, we use a closed formed analytic expression for thickness that can be rapidly computed without the need for computation of streamlines as in [8]. We show in Sect. 2.2, using a 1D analogy, that this closed form expression is equivalent to the result obtained using the streamline method and that the thickness measurement is robust to blurring of the image with a unit-integral kernel. Finally, we present results that compare accuracy and robustness with alternative methods for cortical thickness calculation.

2 Materials and Methods

2.1 The Anisotropic Laplace Equation (ALE) Method

We assume that the brain image has been segmented using a partial classification scheme so that each voxel is assigned a fraction of gray matter (GM), white matter (WM) and cerebrospinal fluid (CSF), with the constraint that the fractions sum to unity [13]. We model cortex as a thin sheet constrained by inner and outer cortical surfaces which are set to temperatures 0 and 1, respectively. We then model the propagation of heat between the two boundary layers of the cortex at equilibrium using the anisotropic form of Laplace’s equation. In our anisotropic model, the diffusion coefficient is assumed to be inversely proportional to the fraction of gray matter in each voxel so that pure white matter and CSF are modeled as perfect conductors. The temperature \(\phi (v,t)\) as a function of spatial location v and time t is given by

$$\begin{aligned} \frac{\partial \phi (v,t)}{\partial t}&=\text {div}\left( \frac{1}{f(v)}\nabla \phi (v,t)\right) ,\,\text {subject to }\phi (v,t)={\left\{ \begin{array}{ll} 0 &{} v\in \partial \varOmega _{inner}\\ 1 &{} v\in \partial \varOmega _{pial} \end{array}\right. } \end{aligned}$$
(1)

where \(\varOmega \) is the domain of the computation, bounded by the closed inner surface \(\partial \varOmega {}_{inner}\) and the closed outer pial surface \(\partial \varOmega {}_{pial}\), and f(v) represents the gray matter fraction at location v. In this formulation it is important that all partial volume voxels that contain cortical gray matter are included within the surfaces that bound \(\varOmega \). The equilibrium solution \(\phi _{\infty }\) for the heat equation with anisotropic flow is given by the anisotropic Laplace equation: \(\text {div}\left( \frac{1}{f(v)}\nabla \phi _{\infty }(v)\right) =0,\) subject to the earlier boundary conditions. Using the calculus of variations, this equation can be reduced to the harmonic energy minimization problem:

$$\begin{aligned} \phi _{\infty }(v)=\underset{\psi }{\arg \min }\int _{\varOmega }\left\| \frac{1}{f(v)}\nabla \psi (v)\right\| ^{2}d\varOmega . \end{aligned}$$
(2)
Fig. 1.
figure 1

(a) Inner and pial surface boundaries overlaid on T1-weighted image; (b) cortical thickness computation based on linked distance overlaid on the gray matter fraction image (LD); (c) the solutions of the Isotropic Laplace Equation (LE) and (d) the proposed ALE method are shown as a color coded temperature distribution with green lines depicting the streamlines. The yellow arrow depicts an example region where the three methods differ significantly.

Full size image

Figure 1 illustrates the color coded temperature distribution solution to the ALE in a 2D section of cortex in comparison to the solution of the isotropic LE and the LD method within \(\varOmega \). We also show corresponding streamlines for the ALE using partial tissue fractions relative to those for the isotropic case. Note that the green streamlines in Fig. 1(d) extend into the white matter due to the non-zero gray matter fraction in this region. Integrals of \(\phi _{\infty }(v)\) over these streamlines from inner to pial surface can be used to compute thickness [14]. We propose an alternative simple analytic expression by defining the cortical thickness T(v) at each point on the mid-cortical surface as:

$$\begin{aligned} T(v)&=f(v)\frac{1}{\left\| \nabla \phi _{\infty }(v)\right\| },\,\text {where}\,(v)\in \partial \varOmega {}_{mid}. \end{aligned}$$
(3)

Here the mid-cortical surface \((\partial \varOmega )_{mid}\) is defined as the level set: \((\partial \varOmega )_{mid}=\left\{ v\in {\varOmega }|\phi _{\infty }(v)=\frac{1}{2}\right\} .\) The analytic expression in Eq. 3 is shown in Sect. 2.2 to be equivalent to the path integral for the 1-dimensional solution to the ALE, not only at \((\partial \varOmega )_{mid}\) but at all points with non-zero gray matter fraction. An intuitive explanation of why this approximation works is as follows. We impose boundary conditions of temperatures 0 and 1 on the inner and pial surfaces, respectively. So for homogeneous gray matter the temperature gradient between them, \(\left\| \nabla \phi _{\infty }(v)\right\| \), will be inversely proportion to thickness. Consequently, calculation of the reciprocal of the gradient at the midcortical surface should produce a good estimate of thickness. For the anisotropic case, we account for the increased flux in partial volume voxels that may lie on the mid-cortical surface by scaling by the gray matter fraction f(v).

2.2 Analysis Using a 1D Model

To illustrate how ALE works, we use a 1D model. In this model we assume pure white matter is the region from \(x= -\infty \) to \(-L\), the cortex (containing gray matter) from \(x=-L\) to L, and pure CSF from \(x=L\) to \(\infty .\) We assume an arbitrary gray matter fraction distribution f(x) on \((-L,L)\). We then blur this distribution with a kernel g(x) with the property \(\intop _{-\infty }^{\infty }g(x)=1\). Following Eq. 3, we obtain the thickness

$$\begin{aligned} T=h(x)\left( \frac{d\phi _{\infty }(x)}{dx}\right) ^{-1}\text {at}\ x:\phi _{\infty }(x)=0.5. \end{aligned}$$
(4)

Here \(h(x)=f(x)\) for the unblurred case and \(h(x)=f(x)\circledast g(x)\) for the blurred case, where \(\circledast \) indicates convolution.

The 1D anisotropic heat equation is \(\frac{\partial \phi (x,t)}{\partial t}=\frac{\partial }{\partial x}\left( \frac{1}{f(x)}\frac{\partial \phi (x,t)}{\partial x}\right) \,\) for the unblurred case, with boundary conditions \(\phi (-L,t)\) = 0 and \(\phi (L\), t) = 1. Solving for \(\phi _{\infty }(x)\) gives \(\phi _{\infty }(x)=\intop \limits _{-L}^{x}f(y)dy/\intop \limits _{-L}^{L}f(y)dy\). Substituting this into Eq. 4 gives the correct value \(T=\intop _{-L}^{L}f(y)dy.\) It is somewhat surprising that the thickness value computed in this manner does not depend on the point x at which it is computed in Eq. 4, although for consistency of definition, we always compute it at the mid point between inner and outer boundary.

Now if we blur the gray matter fraction with a unit-integral kernel, solve the heat equation in steady state, and substitute in Eq. 4, the resulting solution is \(T=\intop _{-\infty }^{\infty }\intop _{-\infty }^{\infty }f(x)g(y-x)\,dx\,dy=\intop _{-\infty }^{\infty }f(x)\left( \int _{-\infty }^{\infty }g(y-x)\,dy\right) \,dx=\intop _{-L}^{L}f(x)\,dx\ \int _{-\infty }^{\infty }g(y)\,dy=\intop _{-L}^{L}f(x)\,dx.\) In other words, the thickness calculations using ALE is unaffected by blurring with a sufficiently narrow kernel with unit integral, provided the blurred gray matter fraction of the cortex lies within the bounds defined by the inner and pial surfaces. It should be noted that in the 3D case, blurring of the image will lead to some smoothing of the computed thickness values along the surface, but this should not add bias to the thickness values.

3 Applications and Results

3.1 Average Cortical Thickness Study

We analyzed 3D structural brain MRI scans of 198 normal right handed subjects (76 male, 122 female, age range: 18–26 years) from the 1000 Functional Connectomes Project http://fcon_1000.projects.nitrc.org. Images were acquired using MPRAGE on a SIEMENS TRIO 3T scanner: TR = 2530 ms, TE = 3.39 ms, slice thickness = 1.33 mm, flip angle = \(7^\circ \), inversion time = 1100 ms, FOV = 256 mm \(\times \) 256 mm, in-plane resolution = 256 \(\times \) 192, 128 slices. We used the BrainSuite software (http://brainsuite.org) [13] to define the partial volume fractions, as well as for extracting inner and pial cortical surfaces. All cortical surfaces were aligned to a common atlas space using BrainSuite in order to compute population based averages. Cortical thicknesses were computed as described in Sect. 2.1 where we used a finite difference method to solve the anisotropic Laplace equation.

Fig. 2.
figure 2

(a) Histology based thickness map from Von Economo [15]; (b)–(d)Average cortical thickness maps of left hemisphere. Lateral (upper) and medial (lower) views from N = 198 adult subjects computed using: (b) Anisotropic Laplace Equation (ALE), (c) Isotropic Laplace Equation (LE), and (d) Linked Distance (LD).

Full size image

We processed each of the 198 subject images using three methods: LD, LE and ALE, and mapped these to the atlas to compute the point-wise average cortical thickness on the surface. A surface based Laplace-Beltrami isotropic smoothing [16] of \(\sim \)10 mm fwhm was applied to the thickness estimates in the original subject surface to compensate for discretization and small misregistration errors. We used a robust mean estimate in which outliers (the 5% most extreme values) were first removed for each vertex on the surface. The maps of average cortical thickness estimated are shown in Fig. 2(b)–(d). For comparison we include in Fig. 2(a) a pseudo-colored version of Von Economo’s map of cortical thickness from [15] which is based on histological measurements. The patterns of thickness variation across the cortex are similar for all three methods, however the range of values is quite different. The cortical thickness estimates found using ALE were more consistent with the Von Economo estimates in the parietal, occipital and temporal lobes, while they were different in the frontal lobe. It should be noted that the demographics of the subjects used in the histological study was different than the imaged population. Von Economo and Koskinas used brains of mentally healthy Caucasian subjects, 30–40 years of age [17] whereas the population in the data we used had an age range of 18–26 and were scanned in Beijing. This difference in the demographics, and especially the younger population in our study, may account for some of the increased thickness in the frontal lobe [18]. LE shows higher thickness than the Von Economo estimates everywhere on the cortex and LD showed even higher values. A similar comparison using FreeSurfer based thickness computation method was presented in [19] where it was shown that FreeSurfer also tended to overestimate the cortical thickness, although the pattern of cortical thickness was similar.

Fig. 3.
figure 3

Effect of scanner differences on cortical thickness estimates for four methods: absolute thickness difference between estimates from the Siemens Trio 3T and Siemens Avanto 1.5T scanners, averaged over 5 subjects. ALE = Anisotropic Laplace Equation, LE = Laplace Equation, LD = Linked Distance, FS = FreeSurfer.

Full size image
Fig. 4.
figure 4

Histogram of the average absolute differences between thickness measures for ALE, LE, LD and FS for 1.5T Siemens Avanto and 3T Siemens Trio scanners.

Full size image

3.2 Test-Retest Reliability Study for Multiple Scanners

The purpose of this study was to analyze the effects of scanner differences on cortical thickness measurements in the same subject. The data for this study consisted of 5 normal subjects scanned on two different scanners. All the scans were acquired within three days at the University of Iowa. The first scan was acquired using a Siemens Trio 3T scanner: slice thickness 1 mm, TR 2530, TE 3.99 ms, inversion time 1100 ms, in-plane resolution 1 \(\text {mm}^{2}\) and flip angle \(10^{\circ }\). The second scan was acquired using a Siemens Avanto 1.5T scanner: slice thickness 1.5 mm, TR 26 ms, TE 7 ms, in-plane resolution 1.066 \(\text {mm}^{2}\) and flip angle \(30^{\circ }\). The BrainSuite surface extraction pipeline was executed for all scans. Thickness estimates were computed using the three methods and coregistered to the atlas brain using BrainSuite. In addition, we also executed the FreeSurfer pipeline (Version 5.3.0) on these subjects. Laplace-Beltrami surface based smoothing \(\sim \)10 mm fwhm was applied and thickness differences were computed. The absolute value of thickness difference corresponding to 3T and 1.5T Siemens scanners, averaged over the five subjects was computed for all four methods as shown in Fig. 3. We also show histograms of these differences in Fig. 4. As with the simulation study, the ALE method shows the smallest absolute difference (mean 0.1575 mm, sd 0.1034 mm) between the 3T vs 1.5T scanners relative to LE (mean 0.1786 mm, sd 0.1163 mm) and LD (mean 0.2480 mm, sd 0.1520 mm), and FreeSurfer (mean 0.1679 mm, sd 0.0821 mm).

4 Discussion and Conclusion

Studies of cortical thickness usually compare differences in thickness in homologous areas between two groups, or changes in thickness over time during maturation, aging or disease progression. For this reason, consistency and robustness of thickness estimates is possibly as important as absolute accuracy. For this reason, we examined not only the average cortical thickness over a relatively larger population, but also the consistency of thickness estimates among subjects scanned in two different scanners. The study of 198 subjects produced average thicknesses with the ALE method that are more in line with those reported in the literature using histological studies than the alternative LE and LD methods. Differences in the frontal lobe using ALE compared to the values in the Von Economo atlas are consistent with age differences in the two different populations studied and reported changes in thickness in early adulthood in frontal cortex [18].

The consistency study in Figs. 3 and 4 shows that there is a consistent bias in all methods in the thickness estimates computed from images from the 3T scanner versus those from the 1.5T. However, the histograms of these differences confirm the reduced sensitivity to differences in resolution of the ALE method relative to FS, LE and LD.

In summary, results presented here indicate that ALE is capable of producing cortical thickness estimates that are largely consistent with those reported from histological measurements, and that these estimates are less sensitive to the effects of imaging in a different scanner for the limited range of conditions over which we have so far studied this method. Further evaluation is required, and as with all thickness estimation methods, though inter-scanner differences may be reduced, our results indicate that it is important that scanner-dependent effects be factored into any subsequent analysis.