Regional gray matter atrophy in relapsing remitting multiple sclerosis: baseline analysis of multi-center data.

Regional gray matter (GM) atrophy in multiple sclerosis (MS) at disease onset and its temporal variation can provide objective information regarding disease evolution. An automated pipeline for estimating atrophy of various GM structures was developed using tensor based morphometry (TBM) and implemented on a multi-center sub-cohort of 1008 relapsing remitting MS (RRMS) patients enrolled in a Phase 3 clinical trial. Four hundred age and gender matched healthy controls were used for comparison. Using the analysis of covariance, atrophy differences between MS patients and healthy controls were assessed on a voxel-by-voxel analysis. Regional GM atrophy was observed in a number of deep GM structures that included thalamus, caudate nucleus, putamen, and cortical GM regions. General linear regression analysis was performed to analyze the effects of age, gender, and scanner field strength, and imaging sequence on the regional atrophy. Correlations between regional GM volumes and expanded disability status scale (EDSS) scores, disease duration (DD), T2 lesion load (T2 LL), T1 lesion load (T1 LL), and normalized cerebrospinal fluid (nCSF) were analyzed using Pearson׳s correlation coefficient. Thalamic atrophy observed in MS patients compared to healthy controls remained consistent within subgroups based on gender and scanner field strength. Weak correlations between thalamic volume and EDSS (r=-0.133; p<0.001) and DD (r=-0.098; p=0.003) were observed. Of all the structures, thalamic volume moderately correlated with T2 LL (r=-0.492; P-value<0.001), T1 LL (r=-0.473; P-value<0.001) and nCSF (r=-0.367; P-value<0.001).

RCT Entities:   Population Interventions Outcomes

1. Introduction

Multiple sclerosis (MS) is the most common demyelinating disorder of the human central nervous system. White matter (WM) and gray matter (GM) atrophy and the presence of demyelinating lesions are hallmarks of MS. (Bermel et al., 2003; Cifelli et al., 2002; Ge et al., 2007; Inglese et al., 2007; Henry et al., 2008; Riccitelli et al., 2012; Sharma et al., 2006; Wylezinska et al., 2003). Based on histopathology, atrophy is thought to reflect demyelination, axonal and neuronal loss that occurs early in the disease, but the precise pathological substrate of cortical atrophy remains unknown (Wegner et al., 2006). Both brain atrophy and lesion load are found to be predictors of long term disability (Di Filippo et al., 2010; Popescu et al., 2013).

Magnetic resonance imaging (MRI) is the most common modality for atrophy quantification and lesion load in MS. A number of studies suggest that regional atrophy is a more sensitive indicator of the disease than global atrophy. Deep GM structures appear to atrophy early on in the disease (Audoin et al., 2010; Bergsland et al., 2012; Charil et al., 2007; Fisher et al., 2008; Fisniku et al., 2008; Mesaros et al., 2008; Pagani et al., 2005; Raz et al., 2010; Sabatini et al., 1996; Sepulcre et al., 2006). Deep GM atrophy is found to correlate with clinical disease markers such as fatigue (Niepel et al., 2006), cognition (Brass et al., 2006; Houtchens et al., 2007; Nocentini et al., 2014), expanded disability status scale (EDSS) score (Bakshi et al., 2001; Prinster et al., 2010; Tao et al., 2009) and MR derived measures such as T2 and T1 lesion loads (Ceccarelli et al., 2009; Tao et al., 2009).

Voxel based morphometry (VBM) is frequently used to assess regional atrophy in MS. Studies using VBM have demonstrated atrophy in various regional GM structures (Ceccarelli et al., 2009; Raz et al., 2010; Bendfeldt et al., 2012). In the VBM analysis, each individual scan is spatially normalized to a common stereotactic space, which is then segmented and smoothed to minimize the image registration errors. Segmented volumes are corrected with their respective Jacobian determinant (JD), which represents the volume change. Optimized VBM was implemented using an unbiased template in the VBM analysis. Prinster et al. (2010) implemented optimized VBM and demonstrated atrophy in various deep and cortical GM in relapsing remitting MS (RRMS). Ceccarelli et al. (2012) compared VBM and optimized VBM, DARTEL (http://www.fil.ion.ucl.ac.uk/spm/software/spm8/) in 26 relapsing-remitting MS (RRMS) patients. They also investigated the effect of lesion in-painting prior to segmentation and their results indicate robust atrophy assessment with DARTEL (Ashburner 2007) when combined with lesion in-painting (Chard et al., 2010).

Tensor based morphometry (TBM) is a technique for investigating disease-related changes in brain structures and is shown to provide methodological improvements over VBM (Hua et al., 2013; Kim et al., 2008). The main advantage of TBM over VBM is that the former does not require tissue segmentation. TBM was implemented to investigate structural changes in schizophrenia (Whitford et al., 2007), traumatic brain injury (Kim et al., 2008), dementia (Brambati et al., 2007, 2009), Alzheimer’s disease (Hua et al., 2008a, 2008b, 2013; Leow et al., 2009), HIV/AIDS (Lepore et al., 2008), and pediatric MS (Aubert-Broche et al., 2011). Both VBM and TBM require templates to assess volume changes. Use of an unbiased template to estimate volume changes is shown to be robust against registration errors and improved statistical power (Joshi et al., 2004; Kim et al., 2008; Lepore et al., 2008). An unbiased template can be generated from a group of MRI image volumes using iterative non-linear registration and averaging.

Tao et al. (2009) were the first to implement TBM to estimate regional atrophy in adult RRMS patients. Their results indicated significant atrophy in major regional GM structures that included thalamus, putamen, and caudate nucleus. In that study, the authors reported relatively strong correlation between thalamic atrophy and EDSS and lesion loads. The major limitations of that study were the relatively modest sample size (88 patients), inclusion of patients from a single center, and relatively long disease duration (DD). These limitations also extend to a number of other published studies and limit the generalization of the results. In the present study, we implemented TBM to quantify regional GM atrophy in RRMS patients with relatively short DD. Approximately 91% of patients had DD less than 3 years. These patients were scanned at different centers on different scanners with varying field strengths and manufacturers as part of multi-center clinical trial (Lublin et al., 2013). Such heterogeneity is inevitable in multi-center trials that require large sample size. Multi-center data is critical for identifying MRI-derived measure in evaluating the treatment efficacy and patient management. In this study the regional GM atrophy in MS patients compared to healthy controls was analyzed and the correlation of atrophy with other MRI and clinical scores was estimated on a large cohort. Regional atrophy was estimated using normalized mean JD of the deformation fields obtained following the nonlinear symmetric diffeomorphic registration with an unbiased template (Tao et al., 2009; Avants et al., 2008). We implemented lesion in-painting prior to registration since it was shown to improve the VBM analysis (Ceccarelli et al., 2012; Datta et al., 2014).

2. Materials and Methods

2.1. Patients

All patients included in this study were recruited in CombiRx, a double-blind randomized clinical trial for evaluating the effects of the combination of two drugs, (interferon beta-1a and glatiramer acetate) compared to treatment with either drug alone in RRMS (Clinical trial identifier NCT00211887). Patients were diagnosed using either Poser or McDonald criteria and had more than one relapse in the previous three years (Lindsey et al., 2012; Table 1).

Table 1

Demographic and clinical (for MS only) data of MS subjects and healthy controls

	MS Subjects (924)	Healthy controls (500)
Age (yrs) (mean ± SD, Range)	37.53 ± 9.65, 18–61	37.68 ± 10.93, 20–61
Females (Males)	674 (250)	298 (102)
EDSS (mean ± SD, Range)	1.94 ± 1.15, 0–6	NA
Disease duration (days) (mean ± SD,Range)	583 ± 1168, 16–13318	NA
T2 lesion load (cc) (mean ± SD,Range)	10.60 ± 11.52, 0.06–76.23	NA
T1 lesion load (cc) (mean ± SD,Range)	1.69 ± 2.21, 0.01–15.37	NA

2.2. MR acquisition

The MRI protocol included 2D fluid attenuated inversion recovery (FLAIR), dual echo turbo spin echo (TSE), and pre- and post-contrast T1-weighted images (0.94 mm × 0.94 mm × 3 mm voxel dimension). The 3D T1-weighted images (0.94 mm × 0.94 mm × 1.5 mm voxel dimension) were acquired either with spoiled gradient recalled echo (SPGR) or magnetization prepared rapid gradient echo (MPRAGE) sequences.

2.3. Healthy controls

To better simulate the heterogeneity of scanners used for the MS cohort, we included 3D T1-weighted images of 400 healthy controls from publicly available databases (http://www.cma.mgh.harvard.edu/ibsr/data.html;http://www.oasis-brains.org; http://ida.loni.ucla.edu) and from our center. Of the 400 healthy control scans included in this study, 292 were scanned at 1.5 T and 108 at 3T. The healthy controls were to match the CombiRx data in terms of age, gender, scanner field strength, and imaging sequence (Tables 1, 2a, 2b).

Table 2

2a: Distribution of females and males across the scanner field strength in MS andhealthy controls
	MS			Healthy controls
	1.5 T	3.0 T	All scanners	1.5 T	3.0 T	All scanners
Females	493	181	674	232	66	298
Males	176	74	250	60	42	102
All genders	669	255	924	292	108	400

2.4. Unbiased template

An unbiased template was created from the 3D T1-weighted images on 14 healthy controls (acquired on our Philips 3T scanner). All the images were skull-stripped using BET2, that is part of the FSL library (http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/; Smith 2002), and were corrected for intensity non-uniformity using N4ITK (Tustison et al., 2010). The buildtemplate module in ANTS (http://www.picsl.upenn.edu/ANTS/) was applied to these 14 image volumes to generate the unbiased template used in this study. Briefly, the initial template was obtained by taking the average of all 14 image volumes. Each image volume was then co-aligned with the initial template and the average was calculated again to obtain a refined template. This process continued for a given number of iterations. The template thus obtained is referred to as the unbiased template. According to the developers of ANTS, 10 images are sufficient to obtain an unbiased template (http://www.picsl.upenn.edu/ANTS/). In our study, we generated the unbiased template using 14 image volumes and 8 iterations (Avants et al., 2008).

2.5. Segmentation of T2 and T1 lesions in MS

Dual echo FSE and FLAIR images were processed to segment brain into GM, WM, CSF, and T2-hyperintense lesions using non-parametric and parametric classifiers (Sajja et al., 2006). Briefly, FLAIR, T1 pre- and post-contrast images were co-aligned with FSE images using rigid body registration (Ashburner 2000) followed by skull-stripping the T2 images using the in-house developed semi-automated algorithm (Datta et al., 2006). T2 lesions, parenchyma, and CSF were classified based on the late-echo FSE and FLAIR images using the non-parametric Parzen classification (Duda et al., 2001). The parenchyma was further classified into GM and WM using the parametric technique, expectation-maximization hidden Markov random field (EM-HMRF) (Zhang et al., 2001). False lesion classifications were minimized using the ratio of first-echo and late-echo FSE images. T2 lesions were delineated using fuzzy connectivity technique (Udupa et al., 1998). T1 hypointense lesions were classified using T2, FLAIR, and T1 pre-contrast images (Datta et al., 2006). Hypointense regions within previously classified T2 lesions were identified using morphological grayscale reconstruction technique. T1 post-contrast images were utilized to prevent enhanced lesion from being classified as T1-hypointense lesion. Segmented lesions were corrected for false classifications by experts as described elsewhere (Datta et al., 2006, 2007; Sajja et al., 2006). The MR measures that include T2 and T1 lesion loads (in cc) and normalized CSF (nCSF = CSF volume/whole brain volume) were calculated for each subject. The whole brain volume that includes GM, WM, CSF, and lesion volumes was calculated following the skull-stripping of T2 images.

The 2D T1 pre-contrast images were co-aligned with the 3D T1 images using nine parameter affine tranrmations. The transformation matrix was applied to segmented images that contained T2 and T1 lesions to deform them to the 3D T1 space using nearest neighbor interpolation. An automated pipeline was developed to obtain the transformation matrix and applied to deform the segmented images of all the patients.

2.6. Lesion in-painting

Lesion in-painting in MS is thought to improve the quality of nonlinear registration (Sdika and Pelletier, 2009; Chard et al., 2010; Ceccarelli et al., 2012). Lesion in-painting was applied to the T1-weighted images following the classification of T2 lesions. Intensities of voxels on the T1-weighted images that correspond to T2 lesions on the FLAIR images were replaced by WM intensities. Briefly, the 2D FSE and the corresponding tissue/lesion segmented images were aligned with 3D T1 image using affine transformation. The full width at half maximum values of global WM intensity in T1 images were obtained using the voxel intensity of WM tissue as classified using PD and T2 images. Intensities of lesion voxels were then randomly replaced with the intensity from a set of global WM intensity values. Gaussian filter was applied to minimize the noise due to random distribution of these intensities. The final in-painted images appear as T1-weighted images with normal WM (Datta et al., 2014).

2.7. Regional atrophy

The T1-weighted images of MS and healthy controls were reformatted to axial orientation followed by skull-stripping using BET2 module of FSL as described above. All the images were co-registered to the unbiased template using non-linear symmetric diffeomorphic registration. An automated pipeline that incorporated these modules was implemented. All images were re-sliced to obtain 1 mm3 isotropic voxels. Parameters associated with BET2 were optimized on 50 scans and these parameters were fixed for all the scans included in the study. Atrophy at each voxel was calculated using the logarithm of the JD, obtained from the 3D diffeomorphic maps that were created during non-linear registration with the template. The logarithmic of JD at each voxel was normalized as described elsewhere (Tao et al., 2009). Positive and negative values of the normalized logarithm of JD represent loss of tissue indicating atrophy and enlargement of tissue indicating hypertrophy, respectively.

The structures associated with the ICBM template (http://www.loni.ucla.edu/ICBM/Downloads/Downloads_ICBMtemplate.shtml) were used to automatically identify regions on the deformed images. First, the ICBM template was co-registered to the unbiased template using symmetric diffeomorphic non-linear registration. The deformation field thus obtained was applied to the structural image of the ICBM template to map the regional structures onto unbiased template using the nearest neighbor interpolation. The ICBM template has 56 labeled structures including CSF and single WM structure. Also 12 structures (reticular, ventral anterior, dorso-medial nuventral, latero-dorsal, lateral geniculate, lateral posterior, ventral lateral, paraventricular, centromedian, medial geniculate, anterior, and ventral posterior nuclei) were combined to obtain thalamus. Therefore the final number of deep and cortical GM structures was 45. The deformed structural image was applied to images for atrophy measurements of various GM structures. Regional atrophy was measured as the mean of normalized logarithm of the JD within a given structure. Regional volumes of each structure in each patient were evaluated from normalized mean Jacobians as described elsewhere (Aubert-Broche et al., 2011).

2.8. Statistical analysis

TBM Analysis

Analysis of covariance (ANCOVA) was used to assess the atrophy of various GM structures in the MS patients relative to the control cohort using SPM8. Specifically, ANCOVA was applied on normalized logarithm of the JDs to analyze the group differences between controls and MS patients. Age and gender were included in the analysis as covariates and corrected for multiple comparisons using family-wise error (p = 0.05) We repeated the analysis with various cluster sizes (0, 10, 20, 30, 40 and 50 voxels) and did not observe any differences(Tao et al., 2009; Narayana et al., 2010). The regional atrophy between MS patients with EDSS of 0 (n = 118) and healthy controls (400) was analyzed to assess brain atrophy in patients without clinical disability. In order to investigate the effect of clinical disability on regional atrophy, we grouped MS patients into different groups and TBM analysis was performed : 1) between MS patients with EDSS ≤ 3.0 (n = 812) and MS patients with EDSS ≥ 3.5 (n = 112); 2) between MS patients with EDSS ≤ 2.5 (n = 725) and MS patients with EDSS ≥ 3.0 (n = 199); and 3) between MS patients with EDSS ≤ 2.0 (n = 597) and MS patients with EDSS ≥ 2.5 (n = 327).

Effects of age, gender, scanner field strength, and imaging sequence

General linear regression analysis was applied to investigate the effect of age, gender, scanner field strength and imaging sequence on volumes of regional structures in MS. Healthy control data was used to correct for these factors in MS.

Correlations with lesion loads and nCSF

Additionally, correlations of EDSS, T1 lesion load (T1 LL), T2 lesion load (T2 LL), and nCSF with estimated regional volumes were analyzed using Pearson’s correlation coefficients.

3. Results

3.1. Demographics data

Of the 1008 patients enrolled in this study, 84 were excluded due to poor image quality, lack of 3D T1-weighted images, and processing errors, including unsatisfactory skull stripping. In-house developed software was used to automatically identify scans with poor quality (Narayana et al., 2013). This resulted in a sample size of 924 at baseline. Figure 1 shows the distribution of age in both MS and age-and gender-matched control groups. The two groups did not differ in age (p = 0.81). No significant difference was observed in the ratio females to males between these two groups (p = 0.60, χ2 test). Demographic data (age, gender, on these patients and healthy controls is summarized in Table 1. Distributions of EDSS, disease duration (DD), T2 LL, and T1 LL of MS patients are also included in the figure. The distribution of gender and imaging sequences across different field strengths are summarized in Tables 2a and 2b.

Figure 1

First row: Age distribution of MS patients and healthy controls, distribution of EDSS in MS; second row: distributions of DD, T2 LL, and T1 LL in MS.

3.2. Non-linear registration

Figure 2 shows examples of non-linear registration of healthy control and MS subject to the unbiased template. Normalized logarithmic JDs are also included in the figures. Lesions in MS subject were in-painted as described earlier. The processing pipeline for analyzing atrophy in MS is summarized in Fig. 3.

Figure 2

Nonlinear registration of healthy control and MS patient to the unbiased template. First row: template, T1 of healthy control, registered T1, and logarithmic Jacobian determinant images; second row: T1 of MS subjects, lesion in-painted T1, registered T1, and logarithmic Jacobian determinant images.

Figure 3

Flowchart of the processing pipeline for atrophy estimation in MS. Here, pre-processing includes skull-stripping and bias correction, lesions comprises of both T2 and T1 lesions with T1 lesions being subset of T2 lesions.

3.3. Regional atrophy

Figure 4 shows the regional atrophy using ANCOVA superimposed on the T1-weighted in MS. The major GM structures that showed significant atrophy are labeled in this figure. Atrophy in the MS patients relative to healthy controls is shown in the first column of Fig. 4. Also included in this figure are the subgroup analyses based on scanner strength and gender. Age, gender, scanner field strength, and imaging sequence were used as nuisance covariates for ANCOVA on groups consisting of whole samples. For comparison of MS and healthy controls scanned at 1.5T and 3T scanners, age and gender were used as nuisance covariates, whereas only age and scanner field strength were considered as nuisance variable for groups consisting of females and males separately. A family-wise error rate of 0.05 was used for correcting for multiple comparisons and a cluster size of 10 was applied on final t-statistic map.

Figure 4

TBM estimates of regional atrophy shown as color blobs superimposed on the unbiased template. The regional atrophy was determined based on the ANCOVA different sub-cohorts. First column: 924 MS vs 400 healthy controls; second column: 669 MS vs 292 controls acquired on 1.5T scanner; third column: 255 MS vs 107 controls acquired on 3T scanner; fourth column: female MS (674) vs female controls (298); and fifth row: male MS (250) vs male controls (102). Age and gender were considered as nuisance variables for ANCOVA analysis in first three groups whereas nuisance variable were age and scanner in analyses between gender groups.

Major structures including thalamus, globus pallidus, putamen, caudate nucleus, brain stem, pons, cerebellum, hippocampus, septal nuclei, insular cortex, and cingulate, temporal, occipital, inferior frontal, parahippocampal, lingual, fusiform, cuneus, and pre-cuneus gyri showed significant atrophy in MS patients compared to healthy controls. This analysis was then repeated separately at 1.5T and 3.0Twith similar results. However, there were few exceptions, superior frontal gyrus, precentral gyrus, superior parietal gyrus, and medulla showed atrophy at 1.5T, but not at 3.0T. Only supramarginal gyrus showed atrophy at 3.0T but not at 1.5T. Notable common structures, thalamus, globus pallidus, putamen, caudate nucleus, brain stem, pons, cerebellum, hippocampus, septal nuclei showed significant atrophy in both female and male MS subjects‥ Figure 5 shows atrophy in 118 MS patients with EDSS score of 0. Significant atrophy was noted in a few structures including thalamus, globus pallidus, putamen, caudate nucleus, brain stem, and cerebellum. No significant difference was observed in the following comparisons: 1) between MS patients with EDSS ≤ 3.0 (n = 812) and ≥ 3.5 (n = 112); 2) between MS patients with EDSS ≤ 2.5 (n = 725) and ≥ 3.0 (n = 199); and 3) between MS patients with EDSS ≤ 2.0 (n = 597) and ≥ 2.5 (n = 327). This suggests that atrophy occurs early on in disease and does not significantly vary with clinical disability as assessed by EDSS.

Figure 5

TBM estimates of regional atrophy as obtained using ANCOVA analysis of MS patients having EDSS 0 (118) compared to healthy controls. Age, gender, scanner field strength, and imaging sequence were considered as nuisance variables for ANCOVA analysis.

3.4. Effects of age, gender, scanner field strength, and imaging sequence andon the volumes of regional structures in MS

Age, gender, scanner field strength and imaging sequence were found to have significant effect on the volumes of 28, 13, 18, and 14 structures respectively in healthy controls. Only 4 structures remained unaffected by these factors. Therefore, the regional volumes of all the structures in MS were corrected for age, gender, scanner field strength and imaging sequence using the regression coefficients from the healthy controls. Generalized linear regression analysis was then applied to the corrected volumes of regional structures in MS. Effects of these factors were estimated by calculating the partial sum of squares, F ratios, and the crresponding p-values. Significant terms with p < 0.001 were identified. Globus pallidus, brain stem, pons, and supramargfinal gyrus were found to be affected by age. Gender had significant effect on thalamus, caudate nucleus, and lingual gyrus. The effect of scanner field strength was found to be significant for caudate nuelues, parahippocampal and supramarginal gyri whereas imaging sequence had effect on brain stem, pons, amygdala, hippocampus, septal nuclei, and postcentrat gyrus.

3.5. Correlations of EDSS, disease duration, T1 LL, T2 LL, and nCSF with regional atrophy in MS

Pearson’s correlations of EDSS, DD, T1 LL, T2 LL, and nCSF with corrected volumes of regional structures are tabulated in Table 3. Only significant correlations (r > 0.1) with p-value < 0.05 were included in the table. Volumes of thalamus, putamen, caudate nucleus, inferior frontal gyrus, anterior commissure, septal nucleus, amygdala, and pons were negatively correlated with T1 LL, T2 LL and nCSF. Moderate correlations with T1 LL, T2 LL and nCSF were observed for thalamic volume. Only hippocampus positively correlated with these measures. Figure 6 shows scatter plots of volumes of thalamus with EDSS, DD, T1 LL, T2 LL, and nCSF. Pearson’s correlation coefficients along with p-values are included in this figure.

Table 3

Correlations of EDSS, disease duration (DD), T2 LL, T1 LL, and nCSF with regional GM atrophy. Only significant correlations with r > 0.1 are included.

Structures	EDSS	DD	T2 LL	T1 LL	nCSF
thalamus	−0.133		−0.492	−0.473	−0.367
putamen			−0.193	−0.176	−0.161
superior frontal gyrus				−0.111	−0.114
cingulate gyrus			−0.123	−0.114
inferior frontal gyrus			−0.120	−0.101
anterior commissure			−0.142	−0.135	−0.210
caudate nucleus			−0.188	−0.173	−0.159
hippocampus		0.120	0.117	0.169	0.130
septal nuclei			−0.146	−0.153	−0.180
subcallosal area					−0.104
mammillary bodies			−0.132	−0.156
red nucleus		−0.118			−0.150
medulla		−0.128			−0.125
amygdala	−0.123		−0.109		−0.127
pons	−0.104	−0.114			−0.201

Figure 6

Scatter plots of volume of thalamus with EDSS, DD, T2 LL, T1 LL, and nCSF in MS patients.

4. Discussion

In this study, we analyzed the regional atrophy of various GM structures in MS in relatively large cohort. We employed TBM, a robust method for assessing the atrophy. Major deep GM structures that include thalamus, putamen, caudate nucleus, globus pallidus, and brain stem and cerebellum showed significant atrophy in MS as compared to healthy controls. The trend is preserved even when the whole group was divided into two groups based on either gender or scanner field strength. The MS group with EDSS 0 showed significant atrophy in these structures but atrophy was not observed within different EDSS groups. These results suggest that atrophy occurs early on in the disease. Generalized linear regression analysis was applied to assess the effects of age, gender, scanner field strength, and imaging sequence on the estimated volumes of regional GM structures. Very few structures were significantly changes with age, gender, scanner field strength, or imaging sequence. Additionally, only thalamic volume was found to be moderately correlated with EDSS, T1 LL, T2 LL, and nCSF. To the best of our knowledge, this is the first study to analyze GM atrophy in a large cohort of RRMS subjects that were recruited in a multi-center clinical study (Bendfeldt et al., 2012). The patients were systematically evaluated for clinical and MR scans using strict protocol.

Our studies indicate that thalamus is one of the few structures that consistently showed atrophy across gender and scanner field strength. This along with previously published studies suggests that thalamic atrophy is a robust measure in MS and could play significant role in evaluating the efficacy of treatment. Based on histopathologic evidence, neuronal loss in thalamus may be responsible for the observed thalamic atrophy (Cifelli et al., 2002; Wylezinska et al., 2003). The supramarginal gyrus showed significant atrophy in males only for which we have no explanation. Overall, atrophy appeared to be more significant in females than males.

As stated above, significant atrophy was observed in thalamus, putamen, caudate nucleus, hippocampus, globus pallidus, red nucleus, cingulate gyrus, and septal nuclei in MS patients even in the absence of any clinical disability (Fig. 5). However, caution should be exercised in interpreting these results because of the relatively small sample size with EDSS of 0. Additionally, we did not observe any significant differences in atrophy with increased EDSS. In, Tao et al. (2009) reported significant atrophy between patients with EDSS ≤ 3.0 and EDSS ≥ 3.5. This difference could be partially attributed to the fact that the disease duration in the present cohort is relatively low.

We analyzed the effect of age, gender, and scanner field strength, and imaging sequence on the estimated regional GM volumes. Our results indicate that these factors have significant effect on the observed atrophy in the majority of the structures. Our results on gender dependence are somewhat different from those reported by Riccitelli et al. (2012) who did not find any gender effect on GM atrophy. Atrophy also appears to be field and sequence dependent in certain structures.

Moderate negative correlation of EDSS with regional volumes of deep GM structures including thalamus, and globus pallidus pars externa were observed. However there was lack of correlation of EDSS with putamen and caudate nucleus. We also observed negative correlation between EDSS and other GM structures including brain stem, entorhinal cortex, medulla, amygdale, and pons. There appears to be some disagreement about the correlation between EDSS and GM atrophy in previously reported studies. For example, significant correlation of EDSS with atrophy in thalamus, putamen and caudate nucleus were reported by Tao et al., (2009). Similarly, Tedeschi et al. (2005) reported significant correlation between GM atrophy and EDSS at the onset of the disease. Shiee et al. (2012) also reported a moderate correlation between EDSS and thalamic atrophy. Our results do not confirm previously reports by Audoin et al. (2010), Ceccarelli et al. (2009), Prinster et al. (2006, 2010), Riccitelli et al. (2012), and Rudick et al. (2009) in which no correlation was observed between regional GM atrophy and EDSS. Additionally, in our studies disease duration showed weak negative correlation with the volumes of thalamus, calcarine fissure, red nucleus, medulla and pons, weak positive correlation with parahippocampal gyrus, hippocampus, and fusiform gyrus.

Both T2 LL and T1 LL correlated negatively with volumes of thalamus, putamen, parahippocampal gyrus, cingulated gyrus, inferior frontal gyrus, anterior commissure, caudate nucleus, amygdala, septal nucleus, mammillary bodies, amygdale and pons, and positively correlated with hippocampus. T2 LL (not T1 LL) correlated with regional volumes of middle frontal gyrus, superior temporal gyrus, fusiform gyrus, and subcallosal area whereas T1 LL (not T2 LL) correlated with regional volumes of superior frontal gyrus, inferior occipital gyrus, lingual gyrus, red nucleus, and medulla. The correlation of thalamic atrophy with lesion load is consistent with those reported by Ceccarelli et al. (2009), Shiee et al. (2012), Prinster et al. (2006), and Tao et al. (2009). In contrast, Riccitelli et al. (2012) did not find any correlation between lesion load and atrophy in thalamus. Normalized CSF (nCSF), a measure of global atrophy was found to correlate with volumes of approximately half of the GM structures however its correlation with thalamic volume appears most significant. Significant correlation between nCSF and thalamic atrophy were also reported previously by Houtchens et al. (2007) and Tao et al. (2009).

Hippocampus showed opposite correlation with T2 LL and nCSF compared to other GM structures (Table 3). While we do not understand the reasons for this opposite behavior, we can speculate that the presence of inflammation and demyelinating lesions in hippocampus that are not detected by conventional sequences (Geurts et al., 2007, Sicotte et al. 2008) may be the reason. Addition of advanced MR sequence such as double inversion recovery might help identify the lesions in hippocampus that could provide better understanding hippocampal atrophy.

The goal of the study was to analyze the feasibility and usefulness of TBM analysis on a relatively large cohort in clinical setting. The processing pipeline including brain extraction, bias correction, and application of ANTS and the creation of normalized logarithmic of the JDs was completely automated. Although the whole procedure was automated, images were checked for errors in brain extraction and mis-registration. This step is the only time consuming step that involves manual intervention. Note that the SPM8 analysis did not require any additional intervention other than data input. This suggests that the methodology, while laborious, can be adapted for analyzing regional atrophy in clinical settings.

The study has several limitations. A major limitation of the study is its cross-sectional nature. In the near future, we plan longitudinal analysis that could help gain better insight into the disease progression and the association between changes in atrophy in various GM structures and clinical measures. Our analysis included EDSS as a measure of clinical disability in this study since the Food and Drug Administration recognizes EDSS only even though the CombiRx data includes Multiple Sclerosis Functional Composite (MSFC) and its components. We plan to analyze the correlation of MSFC and its components with regional atrophy in our future studies. We included 400 healthy controls to assess the atrophy in MS and it is worth noting that these controls were not scanned on the same scanners as MS subjects. Therefore, the effect of different scanner type on the atrophy assessment remains unknown and requires further investigation by acquiring significantly large number of controls and MS subjects on a single and multiple scanners. Scans of healthy controls acquired at 3T scanner were chosen for generating the template. To the best of our knowledge, there is no literature suggesting the effect of field strength on the template. Not all the 3D T1-weighted images were acquired in the axial format, and therefore sagittal and coronal images were reformatted to axial format prior to co-registration with the template. This requires interpolation which introduces some partial volume averaging effect but it is common in image processing. As observed from figure 4, significant differences were observed in WM regions in MS compared to healthy controls. Current study focuses exclusively on GM atrophy. Our future plans include the investigation of WM atrophy and its relation with clinical measures.

In summary, this comprehensive study on regional GM atrophy in a large RRMS data acquired on multiple scanners operating at different field strengths establishes the feasibility of analyzing the data acquired in multi-center study.