Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Med Image Comput Comput Assist Interv. Author manuscript; available in PMC 2017 December 8.
Published in final edited form as:
PMCID: PMC5722235

Gray Matter Surface based Spatial Statistics (GS-BSS) in Diffusion Microstructure


Tract-based spatial statistics (TBSS) has proven to be a popular technique for performing voxel-wise statistical analysis that aims to improve sensitivity and interpretability of analysis of multi-subject diffusion imaging studies in white matter. With the advent of advanced diffusion MRI models – e.g., the neurite orientation dispersion density imaging (NODDI), it is of interest to analyze microstructural changes within gray matter (GM). A recent study has proposed using NODDI in gray matter based spatial statistics (N-GBSS) to perform voxel-wise statistical analysis on GM microstructure. N-GBSS adapts TBSS by skeletonizing the GM and projecting diffusion metrics to a cortical ribbon. In this study, we propose an alternate approach, known as gray matter surface based spatial statistics (GS-BSS), to perform statistical analysis using gray matter surfaces by incorporating established methods of registration techniques of GM surface segmentation on structural images. Diffusion microstructure features from NODDI and GM surfaces are transferred to standard space. All the surfaces are then projected onto a common GM surface non-linearly using diffeomorphic spectral matching on cortical surfaces. Prior post-mortem studies have shown reduced dendritic length in prefrontal cortex region in schizophrenia and bipolar disorder population. To validate the results, statistical tests are compared between GS-BSS and N-GBSS to study the differences between healthy and psychosis population. Significant results confirming the microstructural changes are presented. GS-BSS results show higher sensitivity to group differences between healthy and psychosis population in previously known regions.

Keywords: NODDI, brain microstructure, Gray matter surface based analysis


Gray matter based spatial statistics method (GBSS) was first proposed to study the cortical microstructure development of preterm human brains [1]. This method adapted the tract-based spatial statistics (TBSS) [2] framework to overcome partial volume contamination in diffusion measures within cortical microstructure where spatially transformed individual cortical measurements are projected onto the mean cortical skeleton by searching in a direction perpendicular to the cortical skeleton identifying maximally probable cortical voxels. Unlike TBSS - where projection was based on high fractional anisotropy (FA), GBSS seeks voxels with high probability of GM. This approach was customized as NODDI-GBSS (will be called as N-GBSS in rest of the paper for ease of reference) in a recent study [3] to take advantage of indices provided in NODDI in segmenting the GM instead of using T1 weighted images. GM concentrations are then used along with FA to construct “T1-like” images in diffusion space that are then used to normalize to standard space (Montreal Neurological Institute, MNI). Thinning of the voxels based on TBSS approach is advantageous in improving the sensitivity of diffusion measures in GM; however, it could lead to missing some of the significant regions of interest within GM. Also, the registration approach followed in N-GBSS has the possibility of compromising the sensitivity of voxel-based correspondence when compared to using better registration and segmentation methods. Caveats associated with using TBSS based approach have been presented in recent literature [4] could hold true for GBSS approaches as well, e.g., anatomical inaccuracies and bias in skeleton projection step, effect of varying the threshold in skeletonization process, noise dependency on the shape of the skeleton.

In this work, we present the GM surface based spatial statistics (GS-BSS) approach where the GM central surface is derived from structural scans using Multi-Atlas Cortical Reconstruction Using Implicit Surface Evolution (MaCRUISE) [5]. MaCRUISE has been shown to have superior performance in achieving surface accuracy by combining multi-atlas segmentation with cortical surface reconstruction. Structural images are segmented and normalized to MNI template space with diffeomorphic anatomical registration using exponentiated lie algebra (DARTEL) method [6] and smooth correspondence between cortical surfaces is obtained with diffeomorphic spectral matching [7]. This method retains the speed advantage of spectral matching approach while assuring smooth correspondence between surfaces by performing spectral analysis of the joint graph Laplacian after an initial surface match. To validate our results, we have applied GS-BS on identifying underlying microstructural changes in GM regions between healthy and psychosis groups.

Individuals with psychosis exhibit microstructure changes in prefrontal cortex, including reduced spine density and shorter total dendritic length [8,9]. Recent in vivo studies have demonstrated the ability of NODDI in identifying the group differences between healthy and psychosis groups in white matter [10] and GM regions [3,11]. Using the NODDI model, neurite orientation dispersion index (ODI), a putative marker of dendritic structure and complexity, was calculated and compared between healthy controls and individuals with psychosis.


Data Acquisition and Preprocessing

In this study, we considered a dataset with 114 subjects of whom 47 are healthy and 67 are classified as psychosis group. T1 MPRAGE (256 × 256 mm FOV, 1 × 1 × 1 mm, TE=2ms, TR=8.95 ms and TI=643 ms) and diffusion weighted MRI were acquired on a 3T scanner (Achieva, Philips Medical Systems, Best, The Netherlands) with a 32-channel head coil. The data were acquired at 2.5mm isotropic resolution with FOV of 96 × 96 using DW PGSE-EPI sequence (TR=2.65s, TE=101ms) and Gmax = ~37.5 mT/m. The diffusion scan protocol included a B0 map, two diffusion shells with b-values of 1000 s/mm2 (24 directions), and a HARDI shell with b-value of 2000 s/mm2 (60 directions). For eddy current correction, the diffusion-weighted images (DWI) are first affinely registered to b0 with 12 degrees of freedom using FLIRT in FSL 5.0 [12]. The registration matrix of each DWI is then used to measure patient movement and the gradient table is rotated accordingly.

The flowchart of the approach followed in our method is described in Fig. 1 where registration and segmentation steps are performed on T1 weighted image.

Fig. 1
For each subject, DARTEL registration and MaCRUISE segmentation performed on T1 and all the metrics of interest transformed to MNI space using DARTEL deformations as shown in dotted line. GM surfaces of all subjects in MNI space are then projected onto ...

Each subject’s structural scan was segmented into GM, white matter, and cerebral-spinal fluid (CSF) tissue classes using the VBM8 toolbox ( for SPM12 ( Following segmentation, the GM tissue class images were normalized to a 1.5 mm template image in MNI space (, included with the VBM8 toolbox using the high dimensional DARTEL normalization method [6]. After an initial affine registration of the DARTEL templates to the corresponding tissue probability maps in MNI space (, non-linear warping of the segmented images was then performed to match the corresponding MNI space DARTEL templates (GM, WM). Forward and inverse deformation fields are saved that can be used in further processing.

Individual T1 images are segmented and GM surfaces are derived using MaCRUISE as described in [5] where inner, central and outer cortical surfaces are reconstructed by using the topology-preserving geometric deformable surface model. GM central surface is then normalized to MNI space by applying inverse deformation field to the vertices. As the cortical surfaces derived between subjects do not have predefined correspondence, non-rigid surface registration is required after normalizing to MNI space. The following procedure is followed to account for misalignment of vertices and get voxel-based correspondence needed for group analysis.

Target selection

The Mahalanobis distance (MD) is calculated from each point in a source surface to target surface which is a multi-dimensional generalization of the idea of measuring how many standard deviations away P is from the mean of D. The mean of this metric as shown in Eq. 1 is taken to be the distance measure between surfaces for comparison, where d is the distance metric from surface Y to surface X, μx and Cx are the mean and covariance of X.


Average distances of individual surfaces from all other surfaces are calculated by taking mean of distance measures for that subject. Finally target surface Ts is chosen based on the minimum distance measure where {(x1,yi), (x2,y2),..,(xn,yn)} are the vertices on the target surface.

Surface projection

Let Ss ={(x1,yi),(x2,y2),..,(xk,yk)} be the vertices on each surface S. The Delaunay triangulation, DT (S), based on the convex hull of points in S [13] is generated which ensures that the circumcircle associated with each triangle contains no other point in its interior. N-D nearest point search is performed to obtain initial match K1, between source surface S and the target surface Ts using DT (S). Similarly, corresponding indices K2, from Ts to S are captured. K1 and K2 are used as starting point for execution of diffeomorphic spectral matching [7]. Weighted adjacency matrices are calculated for both surfaces as W1 and W2 based on vertices and edge information from each surface. The single correspondence graph is obtained based on the volumetric entity formed by the two surfaces and with associated links obtained through K1 and K2. The spectral matching algorithm yields the final transformation files C12 and C21 that can be used to register Ts to S and vice versa. All the points in source surface Si that match with corresponding point on Ts based on the correlation map C12 obtained in above step be Pn. Individual measurements from GM surface can then be projected onto the target surface by taking an average of intensity values from all the points in Pn. Surface coordinates are converted to actual volume coordinates by applying the inverse transform of the matrix from image volume in MNI space to the vertices.

NODDI processing

The NODDI model was applied and processed using Accelerated Microstructure Imaging via Convex Optimization (AMICO) [] to estimate intracellular volume fraction (Vic), isotropic volume fraction (Viso) and orientation dispersion index (ODI), which is defined as in Eq. 2 where K is the concentration parameter of the Watson distribution.


ODI is the metric of interest in further analysis for this study. NODDI metrics are first co-registered to intra-subject T1 through the minimally weighted diffusion image using normalized mutual information as a similarity metric. These are then transformed to MNI space by applying forward deformation field obtained from the DARTEL registration step with nearest neighbor interpolation. The nonlinear transform from the spectral matching is applied to the NODDI data to project onto the target surface for group analysis.

N-GBSS processing

GM fraction based on ODI and FA are calculated and a corresponding pseudoT1 is generated in diffusion space used for the registration to MNI space as proposed in N-GBSS [3]. A threshold of 0.55 is empirically set for GM mask relaxing it to capture more number of highly probable GM voxels based on our dataset as original threshold of 0.65 is yielding further thinner skeleton.

Experimental setup/Statistical Analysis

The proposed method is tested on schizophrenic population in identifying the significant regions of decreased ODI between healthy and psychosis groups. In this study, we considered a dataset with 114 subjects of whom 47 are healthy and 67 are classified as psychosis group. To compare the sensitivity retained after registration in identifying significant regions between the groups, VBM analysis is performed in SPM. All individual ODI maps are merged into a 4D volume to prepare for voxel-based analysis within GM surface skeleton regions. A design matrix is formed with unpaired 2 sample t-tests between the groups adjusting for age. Permutation-based statistics are performed on skeletonized ODI for each method using FSL randomize (10,000 permutations) and results are presented for prefrontal cortex regions of anterior cingulate gyrus (ACG) and median frontal cortex (mFC) in Table 1. P-values < 0.05 family-wise error (FWE) corrected for multiple comparisons are considered to be significant.

Table 1
Significant regions of between group analysis with healthy>psychosis using threshold free cluster enhancement with FSL randomize cluster analysis.


N-GBSS based skeleton and GS-BSS based skeleton are overlaid on mean ODI (Fig. 2). As can be noticed from the qualitative image, the coverage of N-GBSS based skeleton is limited especially in the cortical folding regions. Mean voxel count after skeletonization is around ~56K from N-GBSS where as it is around ~101K for GS-BSS. However, N-GBSS captures subcortical regions (e.g., putamen and caudate) which are not overlapped with cortical surface based skeleton from GS-BSS.

Fig. 2
Mean ODI image with an overlay of N-GBSS based skeleton (blue) and GS-BSS based skeleton (yellow) on different axial slices with corresponding z index in MNI space. GS-BSS closely follows the cortical structure with ~101K highly probable GM voxels compared ...

In the ROI analysis, ODI for the medial PFC sub-region is reduced in psychosis patients. To validate the group differences specific to these regions, permutation tests are performed by passing mask of mFC and ACG to both N-GBSS and GS-BSS skeleton images. Significant results from group comparison between healthy and psychosis patients are shown in Fig. 3 and significant values are presented in Table 1. N-GBSS showed significant differences in the anterior cingulate gyrus (ACG) but with fewer voxels in the cluster compared to that of GS-BSS. Also, the medial frontal cortex (mFC) was identified to be significant with p<0.001 using GS-BSS approach whereas none of the voxels showed significance from N-GBSS approach. Fig. 4(a) shows the overlay of mFC region with skeletons from both the methods where it is noticeable that there is little overlap of N-GBSS skeleton on this region that could be causing this issue.

Fig. 3
Significant regions with FWE corrected p<0.05 where ODI in healthy>psychosis. GS-BSS skeleton is represented in green. Results from GS-BSS on row 1 highlighted in red-yellow where ODI in psychosis group is reduced compared to healthy. ...
Fig. 4
Comparison between N-GBSS and GS-BSS methods (a) Segmented GM skeletons with corresponding multi atlas labels. (b) VBM results for whole brain analysis after registration with threshold at 250 voxels (p<0.001 uncorrected) (c) GM skeletons with ...


While N-GBSS yields significant GM voxels that are sensitive to underlying microstructural changes, it still has some limitations. Fig. 4 illustrates issues with possible improvement from GS-BSS. N-GBSS was not able to capture the mFC that has shown to be significant with GS-BSS with high probability (p<0.01). Secondly, GS-BSS yields larger number of highly probable GM voxels on the overall skeleton (Fig 4b) which is indicative of greater sensitivity to capture underlying changes for further analysis. VBM based ODI results are shown in Fig 4c (uncorrected p<0.001 thresholding at 250 voxels).

After applying the correction for multiple comparisons, none of the regions are shown to be significant in N-GBSS where as significant clusters (FWEcorr p<0.005 at cluster level with 1797 voxels) in the prefrontal cortex are seen in the overall ODI with GS-BSS. Lastly, though the results are validated on NODDI derived indices, GS-BSS can be applied to any parameter of interest for analyzing the results in GM.


In this paper, we have presented a novel technique to carry out voxel wise spatial statistics using GM surfaces. We demonstrated the application of this approach to NODDI based metrics in identifying the group differences between healthy and psychosis subjects. We have shown the quantitative differences between the significant regions identified with N-GBSS and GS-BSS. By utilizing the established methods of registration and segmentation, the sensitivity can be retained for performing voxel-wise group analysis compared to N-GBSS method. Moreover, we have shown that using the GM surface as opposed TBSS based skeleton provides more number of highly probable GM voxels that can be helpful in further analysis.

In summary, when compared to N-GBSS, GS-BSS yielded (1) double the number of most probable GM voxels within cortical regions that can be used for further analysis, (2) better sensitivity in VBM analysis with registration based on structural image, and (3) significant associations in regions that are reported to have altered microstructure in psychosis group from post mortem studies. GS-BSS can be applied to any parameter of interest to understand the underlying diffusion microstructural changes in GM.


This work was supported by NIH R01EB017230 & R01MH102266 & Grant UL1 RR024975-01 & Grant 2 UL1 TR000445-06


1. Ball G, Srinivasan L, Aljabar P, Counsell SJ, Durighel G, Hajnal JV, Rutherford MA, Edwards AD. Development of cortical microstructure in the preterm human brain. Proc Natl Acad Sci U S A. 2013;110(23):9541–9546. doi: 10.1073/pnas.1301652110. [PubMed] [Cross Ref]
2. Smith SM, Jenkinson M, Johansen-Berg H, Rueckert D, Nichols TE, Mackay CE, Watkins KE, Ciccarelli O, Cader MZ, Matthews PM, Behrens TE. Tract-based spatial statistics: voxelwise analysis of multi-subject diffusion data. Neuroimage. 2006;31(4):1487–1505. doi: 10.1016/j.neuroimage.2006.02.024. [PubMed] [Cross Ref]
3. Nazeri A, Mulsant BH, Rajji TK, Levesque ML, Pipitone J, Stefanik L, Shahab S, Roostaei T, Wheeler AL, Chavez S, Voineskos AN. Gray Matter Neuritic Microstructure Deficits in Schizophrenia and Bipolar Disorder. Biol Psychiatry. 2016 doi: 10.1016/j.biopsych.2016.12.005. [PubMed] [Cross Ref]
4. Bach M, Laun FB, Leemans A, Tax CM, Biessels GJ, Stieltjes B, Maier-Hein KH. Methodological considerations on tract-based spatial statistics (TBSS) Neuroimage. 2014;100:358–369. doi: 10.1016/j.neuroimage.2014.06.021. [PubMed] [Cross Ref]
5. Huo Y, Plassard AJ, Carass A, Resnick SM, Pham DL, Prince JL, Landman BA. Consistent cortical reconstruction and multi-atlas brain segmentation. Neuroimage. 2016;138:197–210. doi: 10.1016/j.neuroimage.2016.05.030. [PMC free article] [PubMed] [Cross Ref]
6. Ashburner J. A fast diffeomorphic image registration algorithm. Neuroimage. 2007;38(1):95–113. doi: 10.1016/j.neuroimage.2007.07.007. [PubMed] [Cross Ref]
7. Lombaert H, Sporring J, Siddiqi K. Diffeomorphic spectral matching of cortical surfaces. Inf Process Med Imaging. 2013;23:376–389. [PubMed]
8. Garey LJ, Ong WY, Patel TS, Kanani M, Davis A, Mortimer AM, Barnes TR, Hirsch SR. Reduced dendritic spine density on cerebral cortical pyramidal neurons in schizophrenia. J Neurol Neurosurg Psychiatry. 1998;65(4):446–453. [PMC free article] [PubMed]
9. Konopaske GT, Lange N, Coyle JT, Benes FM. Prefrontal cortical dendritic spine pathology in schizophrenia and bipolar disorder. JAMA Psychiatry. 2014;71(12):1323–1331. doi: 10.1001/jamapsychiatry.2014.1582. [PMC free article] [PubMed] [Cross Ref]
10. Rae C, Davies G, Garfinkel S, Gabel M, Dowell N, Cercignani M, Seth A, Greenwood K, Medford N, Critchley H. Deficits in neurite density underlie white matter structure abnormalities in first-episode psychosis. Biological Psychiatry 2017 [PubMed]
11. Nazeri A, Chakravarty MM, Rotenberg DJ, Rajji TK, Rathi Y, Michailovich OV, Voineskos AN. Functional consequences of neurite orientation dispersion and density in humans across the adult lifespan. J Neurosci. 2015;35(4):1753–1762. doi: 10.1523/JNEUROSCI.3979-14.2015. [PMC free article] [PubMed] [Cross Ref]
12. Jenkinson M, Bannister P, Brady M, Smith S. Improved optimization for the robust and accurate linear registration and motion correction of brain images. Neuroimage. 2002;17(2):825–841. [PubMed]
13. Barber CB, Dobkin DP, Huhdanpaa H. The quickhull algorithm for convex hulls. ACM Transactions on Mathematical Software (TOMS) 1996;22(4):469–483.