Integration of multimodal MRI data via PCA to explain language performance.

OBJECTIVE/
METHODS: Neuroimaging research has predominantly focused on exploring how cortical or subcortical brain abnormalities are related to language dysfunction in patients with neurological disease through the use of single modality imaging. Still, limited knowledge exists on how various MRI measures relate to each other and to patients' language performance. In this study, we explored the relationship between measures of regional cortical thickness, gray-white matter contrast (GWMC), white matter diffusivity [mean diffusivity (MD) and fractional anisotropy (FA)] and the relative contributions of these MRI measures to predicting language function across patients with temporal lobe epilepsy (TLE) and healthy controls. T1- and diffusion-weighted MRI data were collected from 56 healthy controls and 52 patients with TLE. By focusing on frontotemporal regions implicated in language function, we reduced each domain of MRI data to its principal component (PC) and quantified the correlations among these PCs and the ability of these PCs to explain the variation in vocabulary, naming and fluency. We followed up our significant findings by assessing the predictive power of the implicated PCs with respect to language impairment in our sample.
RESULTS: We found significant positive associations between PCs representing cortical thickness, GWMC and FA that appeared to be partially mediated by changes in total brain volume. We also found a significant association between reduced FA and increased MD after controlling for confounding factors (e.g., age, field strength, total brain volume). Reduced FA was significantly associated with reductions in visual naming while increased MD was associated with reductions in auditory naming scores, even after controlling for the variability explained by reductions in hippocampal volumes. Inclusion of FA and MD PCs in predictive models of language impairment resulted in significant improvements in sensitivity and specificity of the predictions.
CONCLUSIONS: Quantitative MRI measures from T1 and diffusion-weighted scans are unlikely to represent perfectly orthogonal vectors of disease in individuals with epilepsy. On the contrary, they exhibit highly intercorrelated PCs in their factor structures, which is consistent with an underlying pathological process that affects both the cortical and the subcortical structures simultaneously. In addition to hippocampal volume, the PCs of diffusion weighted measures (FA and MD) increase the sensitivity and specificity for determining naming impairment in patients with TLE. These findings underline the importance of combining multimodal imaging measures to better predict language performance in TLE that could extend to other patients with prominent language impairments.

Introduction

MRI provides a non-invasive method for measuring structural and functional brain changes accompanying various neurological disorders. Advanced MRI techniques increasingly used to explore in-vivo disease biomarkers include T1- and T2-weighted structural images; diffusion weighted images (DWI) and diffusion tensor imaging (DTI); magnetic resonance spectroscopy (MRS); functional MRI (fMRI); susceptibility weighted imaging (SWI) and manganese-enhanced MRI (Obenaus, 2013). Like many others, the field of epilepsy has benefited from the integration of the vast array of MRI techniques available to clinical practice. Within the context of surgical evaluation for refractory epilepsy, structural MRI is now routinely used for the detection of hippocampal volume loss, sclerosis and other pathologies, while DTI is selectively used to delineate the location and direction of white-matter tracts that need to be spared to the extent possible during surgery (Yogarajah and Duncan, 2008). Co-interpretation of MRI data showing widespread structural changes with the broad neurocognitive impairments seen in patients with temporal lobe epilepsy (TLE) has supported the relaxation of the localization-based assumptions in favor of more network-based ones (Saling, 2009).

Neuroimaging studies of language function in epilepsy have relied predominantly on functional data, exploring language lateralization and reorganization (Rodrigo et al., 2008; Hamberger and Cole, 2011; Dym et al., 2011) and prediction of post-surgical language performance (Swanson et al., 2007; Bonelli et al., 2012), including naming and word generation/fluency (Rutten et al., 2002; Sabsevitz et al., 2003; Bonelli et al., 2011). Studies of brain structure have also provided substantial information related to language functioning. Studies based on MRI volumetry have revealed the importance of the hippocampus to language performances, including naming and fluency (Seidenberg et al., 2005). In addition, reduced volumes of the left prefrontal cortex have been associated with impaired fluency in TLE (Hermann et al., 2003). Recently, several investigations have also addressed the role of white matter tracts in language function using DTI data (see Duffau et al., 2008 for a review) and have used fMRI activation patterns to guide DTI analysis of language networks in patients with TLE (Powell et al., 2007). Studies exploring mean diffusivity (MD) and fractional anisotropy (FA) have found associations between the integrity of the uncinate fasciculus and naming (McDonald et al., 2008; Papagno et al., 2011) while those exploring intraoperative mapping have found associations between superior longitudinal fasciculus (SLF) and inferior longitudinal fasciculus (ILF) and phonemic and semantic paraphasias, respectively (Duffau et al., 2002, 2009).

Overall, these approaches have underscored the utility of MRI in delineating which structures contribute to different aspects of language in TLE and in assessing the risk of post-surgical decline based on this anatomical information. However, past literature has often focused on a single modality in isolation, and has rarely attempted to interpret the totality of the MRI data across multiple techniques. In theory, the neurological information captured by each of these modalities is likely to have a great degree of overlap. However, unique information may also be extracted from various MRI modalities that could enhance our ability to understand neurological diseases and cognitive morbidity, including language dysfunction, in epilepsy. Combinant MRI metrics explored to date have shown a promising improvement in clinical utility and predictive ability of MRI data used in conjunction with other modalities (Thesen et al., 2011).

Multivariate approaches have garnered increasing interest from researchers who wish to characterize larger networks impacted by epilepsy. Recent approaches have included graph theoretical analysis, independent component analysis (ICA), and exploratory factor analysis (EFA), all of which rely on extracting information from a correlation or covariance of variables of interest. As reviewed by Bernhardt et al. (2013), graph-theoretical research has successfully characterized the small-world topology (i.e., occurrence of clusters and hubs, short path lengths on average) of various disease networks. Despite the applicability of the covariation-based techniques used across multiple domains (Guye et al., 2010), extant research often focuses on single-domain matrices (e.g., cortical thickness; task-related activation), with very few groups publishing findings on the comparative strength of structural and functional networks across domains (Liao et al., 2013).

ICA's dimension reduction capabilities and ability to identify the principal components of covariance matrices have made it a proven tool for extracting valuable information in multimodal settings. ICA has been successfully used to extract information from complex EEG data (Stern et al., 2009) as well as co-interpreting simultaneously acquired EEG–fMRI data to increase spatiotemporal resolution (Opitz et al., 1999; Ikeda and Toyama, 2000) and source localization (Dien et al., 2003). More recently, researchers have applied ICA to task-based fMRI data to characterize epilepsy networks (Masterton et al., 2013) and investigate changes in language networks (You et al., 2011), reporting significant differences in patients. Research has also shown that differences in scan site, scanner field strength or scanner manufacturer do not obscure ICA and yield consistent results when dealing with principal component vectors (You et al., 2009).

In this study, we quantify the overlap between the information derived from different MRI techniques using PCA and test these components' unique contributions in explaining language performance in a sample of 56 healthy controls and 52 patients with temporal lobe epilepsy. Specifically, we extract the principal components (PCs) that explain the maximal variance of the correlation from the following data domains: i) cortical thickness, gray–white matter contrast (GWMC) and hippocampal volume data obtained from T1-weighted 3-dimensional images and ii) fiber tract FA and MD estimates based on DTI. These MRI variables have all been shown to reveal disease-related pathology in TLE that may affect cognitive functioning (Hermann et al., 2009). We limit our scope to frontotemporal regions in an effort to focus our analyses on regions for which there is strong empirical support of their contribution to language performance. We quantify the associations between MRI-driven PCs and performance in language-based neuropsychological tests using Pearson's correlations. We also use sequential discriminant function analyses (DFAs) to quantify the contribution of DTI measures to the accurate prediction of language impairment in our sample in terms of sensitivity and specificity.

Methods

Participants

This study was approved by the ethical standards committee on human experimentation at University of California, San Diego (UCSD) and completed according to the standards established in the Helsinki Declaration. Written informed consent was obtained from all participants. Fifty-two patients with refractory TLE (ages 19–64) and 56 controls (ages 18–61) were included in the study. Control subjects were recruited through open advertisement and screened for past or present neurological or psychiatric conditions; those with a history of neurological or psychiatric illness were excluded from the study. All TLE patients either were undergoing or had previously completed presurgical evaluation at the UCSD Epilepsy Center and diagnosed with unilateral TLE by a board-certified neurologist (E.S.T. and V.J.I.) with expertise in epileptology. Seizure laterality was identified based on ictal recordings in a video-EEG monitoring unit (with scalp and/or foramen ovale electrodes), seizure semiology, and neuroimaging results. In particular, patients exhibiting dual pathology on MRI (e.g., tumors) were excluded from the sample. Clinical MRI scans were available on all patients and were visually inspected by a board-certified neuroradiologist for detection of unilateral mesial temporal sclerosis (MTS) which was observed in 31 patients (16 patients with left MTS; 15 patients with right MTS; none with bilateral MTS). Patients with contralateral temporal lobe structural abnormalities, lateral temporal seizure focus and structural abnormalities other than MTS were excluded. Demographic information for controls and patients with TLE are presented in Table 1. There were no significant age differences between the groups (p > .05), but controls had significantly higher levels of education compared to patients (p = .003). The absence of significant age differences was particularly important given the possibility of age-related changes in brain structure and function (Steffener et al., 2012, 2013). Both sexes were similarly represented across all three groups χ2(2) = 0.39, p = .82. The percentage of 1.5 vs. 3 T data was also equivalent across groups, χ2(2) = 4.29, p = .12. Mann–Whitney U-tests revealed no significant distribution differences between patients with left and right TLE in terms of continuous clinical variables, i.e., age of seizure onset, seizure frequency and number of medications (all ps > .1). Chi-squared tests showed no significant group difference between left and right TLE patients in terms of the prevalence of positive history of status epilepticus, mesial temporal sclerosis and positive history of febrile seizures (all ps > .2). However, there was a significant reduction in FreeSurfer estimates of hippocampal volumes between groups for the affected (ipsilateral) hemisphere, F(2,105) = 9.01, p < .001 for the right hemisphere and F(2,105) = 37.92, p < .001 for the left hemisphere. Clinical variables for the patients did not exhibit a statistically significant correlation with ipsilateral hippocampal volumes, except for a significant association between longer disease duration and reduced ipsilateral hippocampal volumes, r(52) = .58, p < .001. All patients were right handed.

Table 1

Participant demographics.

Participants	Control	LTLE	RTLE
N	56	25	27
1.5 T	45	17	16
3.0 T	11	8	11
Age	34.8 (11.8)	40.4 (12.4)	38.2 (13.5)
Education	15.0 (1.9)	13.6 (2.2)	13.6 (1.9)
Sex	33 F	14 F	16 F
Age of onset		17.1 (14.0)	16.6 (14.1)
Sz. freq (per mo.)		7.7 (7.9)	4.5 (3.6)
+MTS (on MRI)		16	15

MRI data acquisition and processing

Image acquisition

Magnetic resonance imaging was performed on two General Electric (GE) scanners. The 3 T scans were acquired using a GE Discovery MR750 3.0 Tesla scanner with an in vivo 8-channel phased-array head coil. Image acquisitions included a conventional 3-plane localizer, a T1-weighted 3D customized FSPGR structural sequence (TE = 3.16 ms, TR = 8.08 ms, TI = 600 ms, flip angle = 8°, FOV = 25.6 cm, matrix = 256 × 192, slice thickness = 1.2 mm), and a 30-directional diffusion-weighted sequence (b-value = 1000, TE = 82.9 ms, TR = 8000 ms, flip angle = 90°, FOV = 24.0 cm, matrix = 96 × 96, slice thickness = 2.5 mm, echo-spacing = 588 ms). The 1.5 T scans were acquired using a GE EXCITE HD scanner with an eight-channel phased array head coil. Image acquisitions included a conventional three-plane localizer, GE calibration scan, two T1-weighted customized 3D MPRAGE-equivalent structural scans (TE = 4.8 ms, TR = 10.7 ms, flip angle = 8°, FOV = 25.6 cm, matrix = 256 × 256, slice thickness = 1.0 mm), and diffusion-weighted sequences (single-shot 51-directional echo planar imaging with isotropic 2.5 mm voxels, b-value = 1000, matrix size = 96 × 96, FOV = 24 cm, slice thickness = 2.5 mm). The scans covered the entire cerebrum and brainstem without gaps. Information on the non-linear B0 correction performed on diffusion data from both scanners is provided below. For all scans within the same scanner, the imaging protocol was identical for all participants. For all scans, the patients were seizure-free for a minimum of 24 h.

Image processing

Images in DICOM format were transferred to a Linux workstation for processing. A custom, semi-automated processing stream combining MATLAB and C++ was used. Cerebral parcellations, subcortical volume segmentation, cortical thickness and gray–white matter contrast (GWMC) estimates were accomplished using FreeSurfer v5.1.0 image analysis suite (http://surfer.nmr.mgh.harvard.edu/). Details of the reconstruction, which involves removal of non-brain tissue, Talairach transformation, sub-cortical and cortical segmentation, intensity normalization and topology correction to estimate gray–white matter and gray matter–CSF boundaries, are described in detail in past publications (Dale et al., 1999; Fischl and Dale, 2000; Fischl et al., 1999a, 1999b, 2001, 2002, 2004a, 2004b; Han et al., 2006; Jovicich et al., 2006; Segonne et al., 2004). The cerebral cortex was parcellated into surface ROIs based on sulcal and gyral structure (Desikan et al., 2006). Cortical thickness was defined as the shortest distance from the gray/white boundary to the gray/CSF boundary at each vertex on the tessellated surface (Fischl and Dale, 2000). For GWMC, the cortical matter was sampled at 0.2 mm from the gray–white matter boundary in both directions and gray–white matter contrast was calculated as difference between T1 white matter intensity and gray matter intensity divided by their averaged intensity at each voxel. Consequently, the range for these values was zero to one, with larger numbers representing increased (i.e., better) GWMC. Frontal and temporal lobe ROIs were included in the statistical analysis (see Fig. 1).

Fig. 1

ROIs selected for PCA analysis. Top panels show the cortical ROIs included in the generation of thickness and contrast PCs; bottom panels show the tracts included in the generation of FA and MD PCs. Selected regions lie within the frontal and temporal lobes.

Diffusion weighted images for each patient were first corrected for warping and head-coil related inhomogeneities in signal sensitivity. The images were pre-processed with a five-step procedure that corrected for eddy current, motion and non-linear warping as detailed in Hagler et al. (2009). The steps used for the intensity/distortion correction and registration of the data to the T1 image are provided in detail in Fjell et al. (2012), with the exception that our images were resampled using linear interpolation to 1.875 mm3 isotropic voxels.

Fiber tracts were derived using a probabilistic atlas containing information about the locations and orientations of different fiber tracts as detailed in Hagler et al. (2009). Atlas information was used to estimate probability that a chosen voxel belongs to a particular fiber tract, which became the basis for identifying and labeling thirteen white matter tracts for the entire brain. Fiber ROIs were applied from the atlas to single subjects using both T1-weighted images based on the similarity of diffusion orientations, with an FA threshold of 0.15 (Wakana et al., 2004). For each fiber tract, mean and anisotropic diffusion estimates were calculated using conventional methods and averaged for the entire tract, where FA showed the unidirectionality of the water diffusion [0,1] and MD showed the overall magnitude of diffusion (Fjell et al., 2012). Tracts traversing any portion of the frontal or temporal regions were selected for analysis as tract ROIs (see Supplementary Table S1). These tracts consisted of the following: fornix, superior longitudinal fasciculus, superior corticostriate fasciculus, inferior frontal superior fasciculus of the frontal cortex, striatum abutting the inferior frontal cortex, anterior thalamic radiations, inferior longitudinal fasciculus, inferior fronto-occipital fasciculus, forceps minor, uncinate fasciculus, corpus callosum and cingulum. The tracts can be seen at the bottom panels in Fig. 1, where they are displayed in two groups for easier visualization.

Neuropsychological testing

Vocabulary subtest (Vocab) from the Wechsler Abbreviated Scale of Intelligence — Second Edition (WASI; Wechsler, 1997), Boston Naming Test (BNT; Kaplan et al., 1983), Auditory Naming Test (ANT; Hamberger and Seidel, 2003), and Letter Fluency (LF) and Category Fluency (CF) subtests from Delis–Kaplan Executive Function System (D-KEFS; Delis, Kaplan, and Kramer 2001) were used to quantify language ability across participants. The Vocabulary subtest measures word knowledge and verbal concept formation; BNT and ANT measure the ability to quickly name visually or orally presented stimuli, respectively. LF and CF measure the ability to name words that begin with the same letter or belong to the same semantic category in 1 min. The literature on naming and fluency deficits in TLE even after controlling for education levels is well-established (e.g., Hamberger and Tammy, 1999; Bell et al., 2003) whereas group differences in vocabulary scores are likely to be driven by an interaction of clinical variables and educational achievement. Raw scores from each test were used in analyses and the scores were corrected for age and education by including these variables as covariates.

Statistical analysis

We defined four imaging data domains for the purposes of our analysis: cortical thickness, GWMC, fractional anisotropy and mean diffusivity. We used a surface-based parcellation scheme (Desikan, 2006) to extract averages for cortical thickness and GWMC based on FreeSurfer. We used tractography estimates for fibers that traverse our region of interest and only included fibers that either originate or terminate within a frontotemporal region. We used PCA to extract the main factor with the highest eigenvalue for each imaging data domain. We observed that the factor model explained more than 70% of the variation on average (i.e., communality) for each of the neuroimaging variables. The use of a factor analytic method in a relatively small sample was justified by the overdetermination and high communality observed in our data (Guadagnoli and Velicer, 1988). Moreover, basing our interpretations on the PCs of each data domain minimizes the risk of interpreting unstable patterns related to the idiosyncrasies of the participant sample. To minimize the risk of working with unstable, idiosyncratic factor loadings, we employed stringent factor-loading thresholds adjusted to the observed communality in the data. Consequently, we adopted a 0.7 threshold for all of the PCs, which meets the minimum communality requirements proposed by Guadagnoli and Velicer (1988) for the interpretability of patterns in smaller sample sizes.

ROI-level (thickness and GWMC) and tract-level (MD and FA) correlation matrices in each data domain were analyzed using PCA. Cortical thickness, GWMC, FA and MD were analyzed separately for the frontotemporal regions in each hemisphere. Interhemispheric tracts were excluded from the FA and MD analyses. Given our interest in the PCs, regardless of their relationship to other factors, the unrotated solution was used in each domain. Given the similarity of the factor loading structure, data from both scanning sites (1.5 and 3 T) were used when deriving the PCs for each domain. Correlational analyses were carried out both as zero-order correlations and as first-order correlations after controlling for confounding factors that affect the derived PCs (i.e., age, education, field strength of the scanner). Given multiple comparisons for eight PCs and five neuropsychological measures (i.e., 40 comparisons), a stringent cutoff was used to evaluate statistical significance (p < .001) of brain–behavior correlations, reflecting a 40-fold correction to the family-wise alpha of .05. Significant correlations were followed up with step-wise hierarchical regressions that juxtaposed the variability explained by the PCs with the variability explained by a commonly-referred measure in TLE research, hippocampal volume loss, based on FreeSurfer estimates using the Desikan et al. (2006) atlas.

The predictive power of these PCs was further probed using a sequential (i.e., hierarchical) direct DFA. The goal of DFA is to determine a function that maximizes the separation between groups, hence achieving the greatest possible accuracy in assigning cases to groups. In this study, the functions are linear combinations of demographic and MRI/DTI variables and are called linear discriminant functions (LDFs) (Boslaugh, 2012). Before we set up our analyses, we narrowed our scope to only include the significant brain–behavior correlations (i.e., naming tests and DTI PCs) and confirmed the multivariate normality and homogeneity of variance assumptions in our sample. We dichotomized impairment in naming performance using common clinical cutoffs, where a scaled score less than or equal to 7 was considered impairment for BNT and a raw score lower than 45 was considered impaired for ANT. We then set up three sequential direct DFAs that predicted an individual's BNT impairment (i.e., membership to the “impaired” group). The first DFA included three predictive variables: age, education and total brain volume. In the second DFA, we added hippocampal volume as a predictive variable to the existing three variables to quantify the classification accuracy gained by the inclusion of this commonly used clinical variable. For patients, the volume of the hippocampus in the affected hemisphere was used; for healthy individuals, the average volumes of both hippocampi were used. In the third DFA, we added the PCs for the DTI measures to the existing four variables to quantify the additional accuracy gains resulting from the addition of these measures to the classification model. A single LDF that combined all the predictive variables was generated at each step. We predicted that different combinations of structural (e.g., hippocampal volumes) and DTI data could give rise to different levels of separation between the naming-impaired and healthy groups. Given the absence of a formal criterion and the appropriateness of core DFA metrics in quantifying this notion, p-value, classification accuracy, group distances for the LDF, and the structure loadings for predictor variables for each DFA were used to evaluate the “usefulness” of each DFA and the contribution of each variable to the classification accuracy. Classification accuracy gains resulting from the inclusion of DTI PCs were quantified using the percentage of those correctly classified and LDF distance between groups in each case, where larger metrics indicated better prediction of impairment. Again, given the absence of a formal statistical criterion, any improvement in these metrics accompanied by a statistically significant p-value for the DFA was considered a “useful” improvement in prediction of naming impairment. The sensitivity and specificity of each DFA were then used to construct ROC curves that comparatively illustrate the type I and II errors resulting from the inclusion of DTI data in classification models. This procedure was repeated for the ANT impairment index. Results were evaluated individually and separately for the two field-strength subgroups (i.e., 1.5 and 3.0 T).

Results

Principal component analyses (PCAs)

Dimension reduction using PCA yielded 3 left and 4 right hemisphere cortical thickness factors; 1 left and 2 right hemisphere GWMC factors; 2 left and 3 right hemisphere FA factors; and 2 left and 1 right hemisphere MD factors. In all of the domains, the eigenvalue for the PC was at least 3 times larger than the preceding (secondary) factor. The variance explained by the secondary factors did not exceed 13.9% for any data domain. Statistics for each domain are provided below.

Cortical thickness

PC derived for the cortical thickness of the left hemisphere explained 43.36% of the cortical thickness variation of the frontotemporal region in this hemisphere (KMO = .883). PC derived for the cortical thickness of the right hemisphere explained 40.56% of the cortical thickness variation of the frontotemporal region in this hemisphere (KMO = .882). Items with loadings exceeding 0.7 are provided in Table 2A. Portions of the inferior frontal, superior frontal, middle frontal, lateral orbitofrontal, middle and superior temporal gyri show strong loadings for both hemispheres.

Table 2A

Factor loadings for cortical thickness PCs.

Left hemisphere ROIs	Loading	Right hemisphere ROIs	Loading
Rostral middle frontal	0.8691	Superior frontal	0.8554
Superior frontal	0.8490	Rostral middle frontal	0.8461
Pars triangularis	0.7989	Middle temporal	0.7959
Pars opercularis	0.7910	Superior temporal	0.7897
Superior temporal	0.7876	Caudal middle frontal	0.7869
Caudal middle frontal	0.7800	Pars opercularis	0.7814
Middle temporal	0.7322	Pars triangularis	0.7698
Lateral orbitofrontal	0.7304	Inferior temporal	0.7691
		Lateral orbitofrontal	0.7595

Gray–white matter contrast (GWMC)

PC derived for the GWMC of the left hemisphere explained 77.07% of the contrast variation of the frontotemporal region in this hemisphere (KMO = .969). PC derived for the GWMC of the right hemisphere explained 74.72% of the cortical thickness variation of the frontotemporal region in this hemisphere (KMO = .963). Items with loadings exceeding 0.9 are provided in Table 2B.

Table 2B

Factor loadings for GWMC PCs.

Left hemisphere ROIs	Loading	Right hemisphere ROIs	Loading
Superior frontal	0.9637	Rostral middle frontal	0.9740
Pars opercularis	0.9631	Superior frontal	0.9700
Pars triangularis	0.9608	Pars triangularis	0.9539
Rostral middle frontal	0.9565	Superior temporal	0.9485
Superior temporal	0.9561	Middle temporal	0.9461
Middle temporal	0.9514	Pars opercularis	0.9414
Fusiform	0.9479	Lateral orbitofrontal	0.9387
Inferior temporal	0.9409	Inferior temporal	0.9328
Caudal middle frontal	0.9341	Pars orbitalis	0.9248
Lateral orbitofrontal	0.9275	Caudal middle frontal	0.9142
Pars orbitalis	0.9045	Fusiform	0.9132

Mean diffusivity and fractional anisotropy

PC derived for FA of tracts in the frontotemporal region explained 49.23% of the existing variation (KMO = .826) in the right hemisphere tract FA and 51.98% of the existing variation (KMO = .851) in the left hemisphere tract FA. Items with loadings exceeding 0.7 are provided below in Table 2C.

Table 2C

Fractional anisotropy (FA) PC factor loadings.

Left hemi. tracts	Loading	Right hemi. tracts	Loading
Left parietal SLF	0.8821	Right IFOF	0.8706
Left temporal SLF	0.8744	Right temporal SLF	0.8438
Left IFOF	0.8586	Right parietal SLF	0.8388
Left uncinate	0.7695	Right uncinate	0.7796
Left ILF	0.7648	Right IFSFC	0.7659
Left SIFC	0.7516	Right ILF	0.7585
Left IFSFC	0.7428	Right SCSPC	0.7135
Left SCSPC	0.7159
Left SCSFC	0.7144

PC derived for MD of tracts in the frontotemporal region explained 81.92% of the existing variation (KMO = .914) in the right hemisphere tract MD and 82.71% of the existing variation (KMO = .925) in the left hemisphere tract MD. Items with loadings exceeding 0.7 are provided below in Table 2D.

Table 2D

Mean diffusivity (MD) PC factor loadings.

Left hemi. tracts	Loading	Right hemi. tracts	Loading
Left SCSFC	0.9741	Right SCSFC	0.9718
Left SCSPC	0.9737	Right parietal SLF	0.9693
Left temporal SLF	0.9672	Right temporal SLF	0.9687
Left parietal SLF	0.9661	Right SCSPC	0.9669
Left IFSFC	0.9649	Right IFOF	0.9443
Left uncinate	0.9415	Right CING (cingulate)	0.9434
Left IFOF	0.9309	Right ATR	0.9299
Left SIFC	0.9138	Right uncinate	0.9168
Left CING (cingulate)	0.9077	Right ILF	0.9139
Left ILF	0.9044	Right IFSFC	0.9087
Left ATR	0.8986	Right pyramidal	0.8919
Left pyramidal	0.8951	Right SIFC	0.8678
Left CING (parahippocampal)	0.8667	Right CING (parahippocampal)	0.8415

Differences in language performance

Omnibus tests showed group differences in raw scores across all language tests, driven by deficits observed in the left TLE group. The results, summarized in Table 3, are consistent with past literature and show that the group differences for all but ANT (p < .05) remain statistically significant after correcting for age and education (p < .001).

Table 3

Language performance across groups.

	Vocab	BNT	ANT	LF	CF
Controls	64.4 (8.7)	54.5 (3.3)	48.7 (1.6)	46.4 (12.7)	44.6 (7.3)
LTLE	45.5 (7.3)	43.5 (7.9)	41.7 (6.5)	27.8 (8.7)	28.8 (9.2)
RTLE	59.0 (7.3)	52.0 (6.5)	46.9 (8.9)	35.4 (12.7)	33.2 (10.0)
Max score	80	60.0	50.0	–	–
Significance	**/**	**/**	ns/ns	**/**	**/**

Asterisks before and after the “/” sign indicate significant levels before and after controlling for age and education, respectively, where ns: not significant; **p < .001.

p < 0.001.

Correlational analyses

Correlations among neuroimaging PCs

Multiple significant correlations within and across domains were identified among the eight PCs representing the four neuroimaging data domains (see Supplementary Table S2) even after correcting for multiple comparisons using p < .005. Left and right hemisphere PCs were significantly correlated in all data domains (p < .0001). Zero-order correlations indicated that increased thickness PCs were associated with increased FA PC and that increased contrast was associated with increased FA and MD. However, these associations were no longer significant after controlling for age, education, field strength differences and total brain volume (i.e., confounding factors). A significant correlation was also observed demonstrating that decreased FA PC was associated with increased MD PC only after controlling for the effect of the confounding factors.

Correlations with language performance

Associations between the PCs and the language measures were derived after controlling for age, education, scanner field strength and total brain volume (i.e.) confounding variables. The partial correlations are summarized in Table 4. The results showed a significant positive association between BNT performance and left FA PCs even after controlling for the effects of confounding factors (p < .001). This association was not statistically significant for the right FA PC. Higher left MD PC was associated with worse ANT performance after controlling for the confounding factors (p = .001). This association was not statistically significant for the right MD PC. Higher left FA PC was associated with improved performance on the WASI-Vocabulary tests even after controlling for confounders; however, given the high number comparisons in Table 4 (i.e., 40), associations weaker than p = .001 were not probed further. None of the cortical thickness or contrast PCs showed statistically significant associations with the language performance tests after controlling for confounding factors. In addition, there were no statistically significant correlations between any of the neuroimaging variables and measures of letter or category fluency.

Table 4

Pearson's correlations between PCs and language test performances.

	Vocab	BNT	ANT	LF	CF
N	56	75	40	73	73
LH thickness PC	−.034	−.030	.257	−.054	.006
RH thickness PC	−.061	−.118	.195	−.095	−.012
LH GWMC PC	.044	.217	.022	.193	−.013
RH GWMC PC	.000	.126	−.011	.125	−.017
LH FA PC	.386*	.403**	.337	.188	.239
RH FA PC	.326	.256	.202	.073	.154
LH MD PC	−.326	−.149	−.517**	−.111	−.101
RH MD PC	−.235	−.035	−.425†	−.052	−.028

Correlations after controlling for confounding factors (i.e., age, education, field strength differences and total brain volume).

p < .01.

p < .005.

p < .001.

Significant correlations for BNT and ANT were probed further by follow-up hierarchical linear regressions with all of the confounders entered in the first block, hippocampal volume entered on the second block, and the appropriate PC entered stepwise in the third block. For BNT, PC for FA accounted for an additional 7.8% of the variation in BNT scores, ΔF(1,73) = 9.28, β = .347, p = .003, over and above the 15.2% explained by hippocampal volume. For ANT, PC for MD accounted for an additional 9.9% variation in ANT scores, ΔF(1,38) = 5.964, β = −.389, p = .019, over and above the 21.1% explained by hippocampal volume.

DFA results and ROC curves

Given their significant correlation with left FA PC and left MD PC, BNT and ANT impairment was further analyzed. Three DFAs were performed to predict membership to the BNT-impaired and the ANT-impaired group separately. The null hypothesis for the homogeneity of covariance matrices was confirmed using Box's M, which was insignificant for all DFAs. Results are summarized in Table 5, which shows the significance level, percent accurately classified by each DFA, group means (i.e., centroids) on the LDF and the structure loadings for each predictor variable (i.e., the correlation between the predictors and the LDF) separately for DFAs predicting BNT and ANT impairment. The highest structure loading in each block (e.g., DFA 1) appears in bold; statistically significant loadings appear in italics. Results indicated that BNT prediction could be predicted with 74.6% accuracy by weighing years of education and left hemisphere FA PC for classification (DFA 3) more than any other classifier (i.e., ipsilateral hippocampal volume) in the LDF. Addition of hippocampal volume (DFA 2) to the BNT classification model decreased accuracy by 1.7% compared to DFA 1; addition of DTI PCs increased accuracy by 6.8% compared to DFA 1. On the other hand, ipsilateral hippocampal volume appeared to be the best classifier for ANT impairment regardless of the inclusion of DTI PCs in the model. Nonetheless, improvements in classification accuracy were still observed for DFA 6 that included DTI PCs. Addition of hippocampal volume (DFA 5) to the ANT classification model increased accuracy by 12% compared to DFA 4; addition of DTI PCs increased accuracy by 27.8% compared to DFA 4. Group centroids on the LDF had negative values for the naming-impaired group for each DFA and positive values for the unimpaired group. Directional correspondence between the centroid values and variables loadings showed that the likelihood of naming impairment decreased for all DFAs (1–6) as education, hippocampal volumes and anisotropic diffusion increased. Distances between group centroids for each DFA indicated that naming-impaired and unimpaired subjects were better classified (i.e., more distant group centroids) for ANT than BNT.

Table 5

Comparison of DFAs that predict BNT and ANT impairment.

	BNT impairment (1.5 T; N = 59)			ANT impairment (3.0 T; N = 25)
	DFA 1	DFA 2	DFA 3	DFA 4	DFA 5	DFA 6
Significance	.017	.027	.004	.152	.042	.015
Correctly classified	67.8%	66.1%	74.6%	68.0%	80.0%	95.8%
Group centroids	−.642 ↔ .305	−.671 ↔ .319	−.999 ↔ .475	−1.013 ↔ .253	−1.490 ↔ .372	−2.924 ↔ .585
Age	−.041	−.039	−.026	.039	.026	−.087
Education	.880	.842	.565	.527	.358	.166
Total brain volume	.439	.420	.282	.588	.400	.232
Ipsilateral HC volume		.483	.324		.717	.605
Left hemi FA PC			.565			.289
Right hemi FA PC			.421			.136
Left hemi MD PC			−.143			−.331
Right hemi FA PC			−.046			−.166

BNT results lost their significance but retained their structure when using the smaller subgroup (i.e., 3.0 T) and when both field-strength subgroups were pooled together. ANT results were insignificant for the 1.5 T data. These results are provided in Supplementary Table S3. Cross-validation of results using leave-one-out classification for each observation in our 1.5 T data resulted with 64.4%, 62.7% and 66.1% classification accuracy for DFAs 1–3, respectively. Similar cross-validation for ANT using 3.0 T data resulted with 60.0%, 76.0% and 83.3% classification accuracy for DFAs 4–6, respectively.

ROC curves that visualize the relative classification performance of each DFA are provided below (Fig. 2). Top left quadrant in each figure represents a desirable classification ability with high sensitivity (i.e., low type I errors) and specificity (low type II errors). Both figures demonstrate that the sensitivity and specificity of detecting BNT and ANT impairment can be improved by using PCs calculated from DTI to level above and beyond those obtained by standard models that rely on hippocampal volume loss to predict impairment in TLE.

Fig. 2

ROCs for the prediction on naming impairment. Left panel shows the performance of DFA 1 (solid line), DFA 2 (small dashes) and DFA 3 (short dashes) in predicting visual confrontation impairment within the patients with 1.5 T data. Right panel shows the performance of DFA 4 (solid line), DFA 5 (small dashes) and DFA 6 (short dashes) in predicting auditory naming impairment within the patients with 3.0 T data. Remaining results are provided in Supplementary Table S3.

Discussion

In this study, we integrated information from multiple domains of MRI data to explain language function in patients with refractory TLE. We used four MRI-derived measurement domains (cortical thickness, GWMC, MD and FA) and limited our analysis to the ROI (for thickness and contrast) and tract-level (for MD and FA) characteristics of the frontotemporal regions. We extracted the PCs from our pool of healthy controls and patients with TLE in each domain and quantified the intercorrelations among the PCs and between these PCs and language performance. By pooling data from patients and healthy individuals for structural and diffusion measures for which there is empirical support of group differences (Bernhardt et al., 2010; Kemmotsu et al., 2011), we hypothesized that the PCs would represent disease related change measured by each MRI variable.

Our results confirmed that different MRI data domains convey varying degrees of unique information that are confounded by individual and scanner variables. Past research has shown stronger associations between volume loss and reduced FA in temporal lobe epilepsy (Keller et al. 2012; Scanlon et al., 2013) compared to healthy young adults, who exhibit modest associations between volume and diffusion parameters (Tamnes et al., 2010). In our sample of 108 participants, total brain volumes appeared to be positively associated with both thickness and FA PCs. We found that FA measures show significant associations with cortical thickness and GWMC measures, but that these correlations are partially explained by the differential impact of age, field-strength and total brain volume on these variables. The strong positive association between field strength and GWMC measures appeared especially critical in artificially inflating the FA–GWMC correlation by generating higher GWMC and FA estimates in higher field-strength scanners. On the other hand, the significance of MD–FA correlations was partially attenuated by age, total brain volume and field strength effects. This finding was consistent with earlier reports on the longer-term stabilization of diffusion parameters in patients with TLE following seizures generating negatively correlated MD and FA estimates stemming from the interplay between axonal integrity (e.g., demyelination), gliosis and neuronal loss (Concha et al., 2010).

Despite the intuitive expectation for their occurrence, variations in quantitative gray and white matter measures caused by magnetic field strength differences are still being researched. Chung et al. (2013) recently reported that FA estimates are systematically higher for higher field strength diffusion-weighted acquisitions. In line with Jovicich et al. (2009) and Pfefferbaum et al.'s (2012) reports of local susceptibility, West et al. (2013) recently reported that the overall brain segmentation total brain volume estimates were reliable and identical across both field strengths. However, the magnitude of disagreement between 1.5 and 3.0 T measures was larger for white matter measures than it was for the gray matter, and that deep brain structures, cerebellum and brain stem were especially susceptible. Our findings expand their reports to GWMC measures and comprise an addition to the emerging literature that argues in favor of pooling T1-weighted anatomical data across 1.5 and 3.0 T measures with the use of appropriate correction procedures (Pfefferbaum et al., 2012).

Of the eight PCs from the two domains of MRI data derived from frontotemporal regions, only two of them showed significant correlations with language measures. Lower left hemispheric tract FA was significantly associated with poorer visual naming performance and worse expressive vocabulary. Similarly, higher left hemispheric tract MD was associated with poorer auditory naming performance (Hamberger and Seidel, 2003). In follow-up direct DFAs performed to gauge the marginal explanatory power of the imaging PCs in verbal naming abilities, both FA and MD left hemispheric PCs retained their significance in explaining visual and auditory naming performance, respectively, after controlling for age, education, field strength, total brain volume and hippocampal volume. Inclusion of these DTI PCs not only improved the classification accuracy (Table 5), but also improved the sensitivity and the specificity of the prediction of naming impairment within our sample (Fig. 2). The results were robust to the specification of model parameters: stepwise procedures that involved all 8 PCs yielded similar results that showed only left FA PC and education as significant classifiers for BNT impairment. Thickness and GWMC PCs did not present unique abilities to explain the performance variation observed in these naming tests. In this regard, hippocampal volume, in conjunction with diffusion-based FA and MD PCs, contributed unique information only towards explaining auditory naming deficits in patients with TLE whereas other T1-based thickness and GWMC PCs did not. On the other hand, impairment in visual confrontation naming (BNT) was best classified by heavily relying on education levels and left hemispheric FA (i.e., weighing these metrics heavily on the LDF). In fact, relying solely on ipsilateral hippocampal volume estimates in the absence of DTI PCs appeared to increase both false positives and false negatives when predicting visual naming impairment (DFA 2). This decrease in accuracy contrasted the substantial accuracy gains generated by adding hippocampal volumes to the prediction of auditory naming impairment (DFA 5). DFAs 3 and 6, which combined hippocampal volumes and DTI PCs as classifiers, more accurately predicted both visual and auditory naming impairment. Indeed, following the combination of DTI PCs and hippocampal volume estimates, auditory naming impairment appeared to be predicted with 95.8% accuracy.

We also noted an interesting dissociation between FA and MD PCs in our results that persisted across multiple statistical approaches. Left hemisphere FA appeared specifically more associated with visual confrontation naming than it did with auditory naming, while left hemisphere MD appeared specifically more associated with auditory naming than it did with visual confrontation naming. This dissociation was observed in terms of both partial correlations (i.e., Table 4) and structure loadings (i.e., Table 5) of FA and MD PCs. Given that both tract FA and MD estimates are based on identical diffusion tensors, the divergence of the diffusion results in this manner is counterintuitive. For BNT, we believe that the smaller effect sizes of MD may contribute to the lack of statistical significance in the directionally meaningful correlations observed. Profusion of published FA–behavior relationships in extant neuropsychological literature with concurrent reports of non-significant MD findings (Charlton et al., 2006; Madden et al., 2009; Seghete et al., 2013; Trivedi et al., 2013) further supports this possibility. On the other hand, this theory still does not explain why MD was the only PC we observed to be associated with ANT. Despite recent reports of similar BNT vs. ANT divergence in post-surgical TLE patients that implicate differences in hippocampus-dependence of these tasks (Hamberger et al., 2010), further research is needed to replicate and elucidate this divergence in TLE.

Past research has established the sensitivity of DTI and tractography in probing structural pathology including demyelination, glial hypertrophy, scarring and apoptosis (Obenaus, 2013). Research has already shown that both hippocampal and extra-hippocampal white-matter fiber architecture is affected by refractory TLE (Hermann et al., 2003; Rodrigo et al., 2007; Diehl et al., 2008; Dulay et al. 2009; McDonald et al., 2008; Yogarajah and Duncan, 2008; Keller et al., 2012; Kucukboyaci et al., 2012). Prior estimates for significant correlations between tract specific white-matter FA and language measures, which range from r = .535 for right temporal FA and naming and vocabulary based language performance in children with focal epilepsies (Widjaja et al., 2013) to r = [.65–.80] for right and left arcuate and naming performance (McDonald et al., 2008b), appear directionally consistent with the associations we report. With regard to the magnitudes, we observe a significant attenuation in the diffusion MRI data–language performance relationship after controlling for age, education and hippocampal volumes. Our reliance on partial correlations controlled for these confounding variables may explain the smaller magnitude of our results. Still, all of these findings suggest that DTI can be especially important in augmenting T1-weighted data and that the diffusion metrics obtained may provide additional clinical utility for understanding language impairments at both the regional and the global level.

In summary, we used PCA to simultaneously integrate and interpret MRI data across different modalities. Based on PCs that represent disease-driven change in our participant pool, we report significant correlations between disease-driven cortical thickness and GWMC changes across hemispheres along with associations between FA reductions and MD increases that accompany the refractory disease process. Furthermore, we show that these changes in the diffusion patterns are associated with language performance, where FA reductions are associated with worse visual confrontation naming and MD increases are associated with worse auditory naming. We also find that the FA and MC patterns are highly sensitive and specific in predicting language impairment at the PC level. Our findings highlight the importance of frontotemporal white matter integrity in language performances in TLE and suggest that microstructural measures derived from DTI may be more powerful estimates of cognitive dysfunction than cortical measures derived from T1-weighted volumes. Moreover, our multivariate methods retain the statistical power that is often lost in ROI-based group comparisons.

Several limitations to our study should be noted. First, we utilize a dimension-reduction technique (i.e., PCA) often used with larger samples (e.g., >250). We rely on the high communality we observed within the frontotemporal regions of interest to boost the validity of our results and ensure good recovery of population factors in samples regardless of sample size, level of overdetermination, or the presence of model error (MacCallum et al., 2001). We also work with data that show a high degree of overdetermination (i.e., high number of items per factor), and limit our scope to items with high factor loadings (>0.8), which maximizes the sample-to-population fit for the extracted factors (Velicer and Fava, 1998). Therefore, we believe that the PCs we derive are interpretable independent of the sample size used (Guadagnoli and Velicer, 1988). Working with PCs also allows us to pool data across field strengths and maintain statistical power with minimal problems arising from multiple comparisons. Nonetheless, our methods also have their shortcomings. While we choose to work with the PC of each data domain to capitalize on the stable, reliable portions of the communalities observed, smaller components that were extracted may also have clinical or research significance. Although the eigenvalues we observed strongly suggested a unifactorial structure in each of our MRI domains, it is possible to extract other, smaller factors. These smaller components, which may include overrepresented correlations with hippocampal morphometry and diffusion, may also capture disease-related changes in MRI measurements (see Supplementary Table S4). However, our sample size reduces our ability to accurately recover population characteristics in our sample for smaller components. Moreover, our PCA has relied on the unrotated factor structure in each domain since extracting multiple, maximally orthogonal factors had not been our goal. However, rotated factor structures could be of interest to other researchers who can use larger samples with the aim of extracting multiple, reliable factors from each domain.