Materials and Methods
Study Participants
From our surgical epilepsy imaging database, we selected 15 consecutive patients undergoing presurgical evaluation for drug-resistant focal epilepsy with radiologically apparent (MRI positive) or histologically proven (MRI positive or MRI negative) focal cortical dysplasia who had available 3T MRI structural imaging carried out at the NIH Clinical Center with our standard epilepsy imaging protocol between 2014 and 2019. Patients scanned on a different scanner or with a different imaging protocol were excluded, as were subjects with low image quality or other image artifacts on visual inspection. Data from a control group of 30 healthy volunteers with no previous history of neurologic, psychiatric or other significant medical illnesses that may affect the central nervous system were also included. Data were collected at the Clinical Center of the National Institutes of Health (NIH; Bethesda, MD). All participants were enrolled in a research protocol approved by the Institutional Review Board; informed consent was obtained from all participants.
MR Acquisition Protocol
All participants were scanned with a Philips Achieva 3T scanner in the NIH Clinical Center Radiology Department. Whole brain sagittal 3D MRI scans were acquired as follows: 1) 3D T1 weighted MPRAGE (T1): TR = 6.8–7.2, TE = 3.2 ms, TI = 900ms, flip angle = 90, voxel size = 0.75 x 0.75 x 0.8, acceleration factor 2 in slice direction, acquisition time = 7:02 min; 2) 3D T2 weighted FSE (T2): TR = 2500, TE = 225–245, voxel size = 1 x 1 x 2 or 1 x 1 x 1, acceleration factor 2 in slice and phase directions, acquisition time = 5:03 min; 3) 3D FLAIR: TR = 4800, TE = 271–415, TI = 1600, voxel size = 0.9 x 0.9 x 1, acceleration factor 2 in slice direction, 2.6 in phase direction, acquisition time = 6:10 min.
Preprocessing
All images were initially inspected for excess motion and overall quality. For each individual subject, T
2 and FLAIR images were co-registered to the T
1 image with an affine transformation and a normalized mutual information cost function using the publicly available software package Analysis of Functional NeuroImages (
AFNI)
\cite{Cox1996a}. The resulting images were sampled to the grid of the subject's T
1 image. Registered images were visually assessed for alignment.
Cortical Reconstruction
Cortical reconstruction was performed using standard processes with
FreeSurfer v6.0.0 \cite{Dale1999,Fischl1999}. The T
1 and T
2 weighted images were used as input to FreeSurfer's processing pipeline. Briefly, this processing consists of removal of non-brain tissue using a hybrid watershed/surface deformation procedure
\cite{Ségonne2004}, automated Talairach transformation, segmentation of the subcortical white matter and deep gray matter volumetric structures
\cite{Fischl2004,Fischl2002}, intensity normalization
\cite{Sled1998}, tessellation of the gray matter white matter boundary, automated topology correction
\cite{Fischl2001}, and surface deformation following intensity gradients to optimally place the gray/white and gray/cerebrospinal fluid borders at the location where the greatest shift in intensity defines the transition to the other tissue class
\cite{dale:99,fischl:99}. Once the cortical models are complete, a number of deformable procedures are performed including surface inflation
\cite{fischl:99}, registration to a spherical atlas which is based on individual cortical folding patterns to match cortical geometry across subjects \citep{HBM:HBM10}, and parcellation of the cerebral cortex into units with respect to gyral and sulcal structure
\cite{Desikan2006,Fischl2004}. Results were visually inspected and manually corrected if needed. FreeSurfer output was converted to AFNI compatible format, and cortical surfaces were deformed to a standard mesh consisting of 36,002 vertices per hemisphere.
Lesion Labels
For the subset of patients with MRI positive focal cortical dysplasias, lesions were traced in the volume using AFNI by an experienced neurologist using the T1 weighted image when possible, informed by the T2 and FLAIR images when necessary. For MRI negative patients, the postoperative T1 weighted image was registered to the preoperative T1 weighted image in the same manner as described previously for T2 and FLAIR contrasts, and the resected region was traced using AFNI. The resulting lesion masks were mapped onto the surface at the gray-white junction in the same manner as described when originally sampling the features onto the surface.
Intensity Correction
Intensity correction was implemented using an in-house procedure that performs an edge-aware local normalization, division by a smooth estimate of the local signal mean within each 3D image for each subject. Further mathematical details and examples are provided in the supplementary materials materials
(link to supp). The algorithm were implemented in Python and the source code publicly available on GitHub as part of the
JEM pac
Feature Generation
For each voxel, in addition to the corrected intensity in each MR contrast, we wished to represent the characteristics of the local neighborhood surrounding each voxel in a manner that was insensitive to the local orientation of the cortical sheet. We therefore created a set of rotationally invariant 3D multi-scale local image filters analogous to the multi-scale oriented filters that have been shown to be: 1) an efficient representation of the local statistics of natural image statistics, 2) an accurate model of neurons in the mammalian visual system, and 3) commonly used in computer vision applications to create a robust, efficient, generic basis set for applications such as texture synthesis \cite{Simoncelli2001}. Our rotationally invariant local features consisted of:
- Gaussian: the image blurred by Gaussian kernels of widths 1, 2, and 3 mm;
- Magnitude of the gradient of Gaussians: the magnitude of the first spatial derivative of the image blurred by Gaussian kernels of widths 1, 2, and 3 mm;
- Laplacian of Gaussian: the trace of the spatial second derivative of the image blurred by Gaussian kernels of widths 1, 2, and 3 mm;
- Frobenius Norm of the Hessian: the matrix norm of the spatial second derivative of the image blurred by Gaussian kernels of widths 1, 2, and 3 mm
This resulted in a final 39 feature vector at each voxel consisting of the corrected T
1, T
2, and FLAIR image intensities, as well as the four local features at three spatial scales for each MR contrast. Because well-described features of FCDs have been found at or in close proximity to the gray-white junction, features were then sampled onto the FreeSurfer generated gray-white junction surface. Further mathematical details and example images are shown in the supplementary materials
(features figure). The source code is also publicly available on GitHub as part of the same Python package (
JEM).
Scaling and Dimensionality Reduction
Each resulting surface feature was standardized across the whole cortex within each subject. Dimensionality reduction and whitening was then performed using principal component analysis (PCA) using the
scikit-learn machine learning library in Python
\cite{scikit-learna}. The PCA model was fit to all cortical vertices in all healthy volunteers and applied to the resulting feature set for each subject. Fourteen components were kept (90% explained variance).
Normative Model
We then sought to create a normative model that would characterize the joint probability density function (PDF) of our features across the entire cortex in our healthy volunteer population, allowing for characterization of the expected range of normal cortical variability. Inspection of the histograms and joint scatter plots of the PCA components, although they were uncorrelated with zero mean and unit variance, demonstrated a fair amount of residual structure in the form of variable degrees of skew and kurtosis (overview fig). To allow for more straightforward estimation of outlierness, as well as similarities and differences between cortical vertices, we applied an invertible non-linear transformation to map the features into a latent representation where the PDF follows an approximately multivariate normal distribution using the Rotation-based Iterative Gaussianization (RBIG) procedure \citep{LapaEtal2011}. RBIG is a fast iterative procedure consisting of a sequence of pairs of transformations: 1) a non-linear transformation applied to each of the columns (marginals) of the data matrix, and 2) a linear transformation applied to the entire data matrix. The non-linear column-wise operation is a univariate Gaussianization that converts percentile scores computed using the rank transformation to standard scores (scikit-learn QuantileTransformer). The orthogonal transformation is performed using PCA, with all components retained after each iteration (figure \ref{fig:gaussianize}). Ten iterations were performed. Following this procedure, the histograms and joint scatter plots of each of the resulting components now had an approximately Gaussian distribution, as expected (overview fig).
Although similar to features used in other imaging applications REFSS, this set of features has not been used to our knowledge to describe normal cortical variability or to detect cortical abnormalities. We therefore wished to determine whether our feature set could be used to predict the values obtained using some more commonly used metrics such as surface-based measurements of curvature, sulcal depth, cortical thickness, and gray/white contrast (as calculated using FreeSurfer), or volumetric measures of myelination as in \cite{Glasser_2011}, calculated by dividing the T1 intensity by the T2 intensity sampled onto the surface. Using our feature set as the input, linear regression models (ordinary least-squares) were estimated for each target metric using scikit-learn. Model training and testing was performed in the group of healthy volunteers. Each model was evaluated using a 10-fold cross validation procedure. Each fold was trained on all cortical vertices from 25 randomly selected healthy volunteers and tested on all cortical vertices from the remaining 5 subjects. Performance was evaluated for each model using the coefficient of determination \(r{^2}\) (regression figure goes here), and effect size was reported as in \cite{cohen1988statistical}: \(r=0.1\) as small, \(r=0.3\) as medium, and \(r=0.5\) as large.