Changes of deep gray matter magnetic susceptibility over 2 years in multiple sclerosis and healthy control brain.

In multiple sclerosis, pathological changes of both tissue iron and myelin occur, yet these factors have not been characterized in a longitudinal fashion using the novel iron- and myelin-sensitive quantitative susceptibility mapping (QSM) MRI technique. We investigated disease-relevant tissue changes associated with myelin loss and iron accumulation in multiple sclerosis deep gray matter (DGM) over two years. One-hundred twenty (120) multiple sclerosis patients and 40 age- and sex-matched healthy controls were included in this prospective study. Written informed consent and local IRB approval were obtained from all participants. Clinical testing and QSM were performed both at baseline and at follow-up. Brain magnetic susceptibility was measured in major DGM structures. Temporal (baseline vs. follow-up) and cross-sectional (multiple sclerosis vs. controls) differences were studied using mixed factorial ANOVA analysis and appropriate t-tests. At either time-point, multiple sclerosis patients had significantly higher susceptibility in the caudate and globus pallidus and lower susceptibility in the thalamus. Over two years, susceptibility increased significantly in the caudate of both controls and multiple sclerosis patients. Inverse thalamic findings among MS patients suggest a multi-phase pathology explained by simultaneous myelin loss and/or iron accumulation followed by iron depletion and/or calcium deposition at later stages.

Introduction

Iron is an essential co-factor in many biochemical processes of the brain, including myelin synthesis (Rouault, 2013). Multiple sclerosis is a demyelinating neurodegenerative disorder where brain iron levels are disturbed (Zecca et al., 2004, Hagemeier et al., 2012a), especially in the deep gray matter (DGM) (Stankiewicz et al., 2014), around plaques (Craelius et al., 1982), and in the normal appearing white matter (Paling et al., 2012, Hametner et al., 2013). While iron around plaques has been found predominantly in non-phagocytosing, M1-polarized macrophages/microglia (Mehta et al., 2013), the mechanisms of iron concentration changes in DGM as well as the involved cell types are currently unknown. It has been suggested that iron overload in DGM induces or amplifies neurodegeneration in multiple sclerosis through the enhanced generation of reactive oxygen species (ROS).

Several studies have exploited the exceptional iron-sensitivity of MRI to evaluate tissue iron levels in the DGM of multiple sclerosis patients in vivo. MRI measures such as T2 hypo-intensity, R2* relaxometry, and susceptibility-weighted phase imaging have been found to be correlated with structural damage (Bermel et al., 2005, Zivadinov et al., 2012), advancing multiple sclerosis disease course (Hagemeier et al., 2012b, Blazejewska et al., 2015), and with physical and cognitive measures of disability (Brass et al., 2006, Burgetova et al., 2010, Hagemeier et al., 2013b, Modica et al., 2015). However, most published studies are of cross-sectional design, which provides limited insights on disease-related temporal tissue changes. Although they may provide information about mean group trends, cross-sectional evidence is unsuitable for understanding change, individual differences, and potential mediators of change. Longitudinal studies have several inherent advantages, especially when measuring age-dependent quantities, such as brain iron concentration (Hagemeier et al., 2013a, Ghadery et al., 2015), as each study participant acts as its own control. In particular, longitudinal MRI investigations allow studying time-dependent effects of iron within and between study groups. However, only recently have longitudinal studies been carried out with iron sensitive imaging techniques in multiple sclerosis (Walsh et al., 2014, Khalil et al., 2015, Uddin et al., 2016). The longest follow-up time in these studies was three years (Khalil et al., 2015) and employed R2* relaxometry. A longitudinal study with quantitative susceptibility mapping (QSM) has not yet been published.

In the present longitudinal study we investigated disease-related tissue changes in multiple sclerosis over two years using QSM (Reichenbach et al., 2015). This recently introduced advanced MR imaging technique allows the measurement of subtle changes of the magnetic susceptibility of tissue and is currently regarded as one of the most sensitive and specific techniques for studying tissue iron in vivo in multiple sclerosis (Langkammer et al., 2013, Stuber et al., 2015). However, similar to R2* mapping, tissue susceptibility is also affected by myelin and calcium (Schweser et al., 2011), but their contributions to the voxel susceptibility are opposed to that of iron (paramagnetic iron increases susceptibility; diamagnetic myelin and calcium decrease it), whereas all these substances increase R2*. Hence, the measurement of magnetic susceptibility in multiple sclerosis may provide information on disease-related tissue changes which are complementary to previously reported findings using R2*. Based on the assumption that iron is accumulated as a secondary effect during the course of the disease, we hypothesized that multiple sclerosis patients have higher susceptibility and lower volumes of the DGM than healthy controls, and that this gap is augmented over time.

Materials and methods

Subjects

One hundred and sixty (160) participants were included, consisting of 120 multiple sclerosis patients (relapsing-remitting [RRMS]: 98, secondary-progressive [SPMS]: 22) and 40 age- and sex-matched healthy controls (female-to-male ratio 3:1). Inclusion criteria for multiple sclerosis patients were: RR or SP disease course, age 18–65 years, and Expanded Disability Status Scale (EDSS) between 0 and 6.5. Patients with a relapse and/or steroid treatment within 30 days prior to MRI were excluded. Participants were excluded if they were pregnant or in case of pre-existing medical conditions known to be associated with brain pathology (e.g., cerebrovascular disease, positive history of alcohol dependence). Multiple sclerosis patients were diagnosed using the revised McDonald criteria (Polman et al., 2011) and clinical disease severity was measured using the EDSS (Kurtzke, 1983) (baseline median: 2.5, inter-quartile range: 2.5; follow-up median: 3.0, inter-quartile range: 2.5). Healthy controls were recruited from volunteers with a normal neurological examination and no history of neurologic disorders or chronic psychiatric disorders. The study was approved by the local Ethical Standards Committee at the University at Buffalo and written informed consent was obtained from all participants.

Image acquisition

Participants were imaged at baseline and follow-up with the same 3T GE Signa Excite HD 12.0 scanner (General Electric, Milwaukee, WI, USA) using an eight-channel head-and-neck coil. No major hard- or software upgrades to the MRI system occurred during the duration of the study. Data for QSM were acquired using an un-accelerated 3D single-echo spoiled gradient recalled echo (GRE) sequence with first-order flow compensation in read and slice directions, a matrix of 512 × 192 × 64 and a nominal resolution of 0.5 × 1 × 2 mm3 (FOV = 256 × 192 × 128 mm3), flip angle = 12°, TE/TR = 22 ms/40 ms, bandwidth = 13.89 kHz (Hagemeier et al., 2012b), and a total measurement time of 8:46 min:s. Raw k-space data of each channel were saved for off-line reconstruction of the images. In addition, the following sequences were acquired: dual fast spin-echo (FSE) proton density- and T2 weighted imaging (TE1/TE2/TR = 9/98/5300 ms; echo-train length = 14); FLAIR (TE/TI/TR = 120/2100/8500 ms; flip angle = 90o; echo-train length = 24); spin-echo T1-weighted (T1w) imaging (TE/TR = 16/600 ms); and a magnetization-prepared 3D high-resolution T1w fast spoiled GRE sequence (TE/TI/TR = 2.8/900/5.9 ms, flip angle = 10o, isotropic 1 mm resolution). The 2D sequences were acquired using a 256 × 192 matrix (frequency × phase) and 256 mm × 192 mm FOV, resulting in a nominal in-plane resolution of 1 mm × 1 mm, with 48 gap-less 3 mm slices for whole-brain coverage.

Image analysis

QSM

The processing for QSM is schematically illustrated in Fig. 1. Magnitude and phase GRE images were reconstructed offline on a 512 × 512 × 64 spatial matrix using sum-of-squares and scalar phase matching (Hammond et al., 2008), respectively. To achieve isotropic in-plane resolution, the k-space was zero-padded in phase-encode direction prior to the processing. In-plane distortions due to imaging gradient non-linearity were compensated (Polak et al., 2015). Phase images were unwrapped with a best-path algorithm (Abdul-Rahman et al., 2007), background-field corrected with V-SHARP (Schweser et al., 2011, Wu et al., 2012) (radius 5 mm; TSVD threshold 0.05), and converted to magnetic susceptibility maps using the HEIDI algorithm (Schweser et al., 2012). Magnetic susceptibility was referenced (0 ppb) to the average susceptibility of the brain. This was done under the assumption that a larger reference region would limit the additional inter-subject variability, compared to a smaller reference region. All data processing was performed with in-house developed algorithms in MATLAB (The MathWorks, Natick, MA).

Fig. 1

Schematic illustration of the QSM processing. Images shown correspond to a 41 years old representative healthy control.

DGM segmentation and volumetry

Anatomical DGM regions were segmented with FSL FIRST on 3D T1w images corrected for T1-hypointensity misclassification (Patenaude et al., 2011), using an in-house dilation-based lesion inpainting algorithm. Mean values of susceptibility were calculated in thalamus, caudate, putamen, and globus pallidus. The hippocampus, which was also segmented by FSL FIRST, was excluded from the analysis because the region generally does not show a strong contrast on susceptibility maps and is located close to the brain's surface, which makes it susceptible to processing-related artifacts.

Whole-brain volumes, normalized for subject head size, were measured with FMRIB's SIENAX cross-sectional software tool (version 2.6) at baseline. DGM volumes were normalized using the SIENAX-derived scaling factor. SIENA (Smith et al., 2001) and SIENAX-MTP (Dwyer et al., 2014) were used to perform longitudinal volume assessments of the whole brain, gray matter, white matter, and lateral ventricles. T2-lesions were identified on T2w/FLAIR images, while T1- and contrast-enhancing lesions were depicted using pre- and post-contrast spin echo T1w images. Lesions were segmented using a semi-automated edge detection contouring and thresholding technique, as described previously (Zivadinov et al., 2001).

Statistical analysis

Statistical analysis was conducted using SPSS (23.0, IBM, Armonk, NY) and R (3.2.1, R Foundation for Statistical Computing, Vienna, Austria). Outcome measures were averaged between both hemispheres and tested for normality using Q-Q plots and the Shapiro-Wilk test. Baseline characteristics were compared using chi-squared and Student's t-test. A mixed factorial ANOVA was performed as the primary analysis. In these models, baseline and follow-up MRI measures (susceptibility and normalized volume) were compared and tested for (i) the effects of time [baseline vs. follow-up], (ii) study group [controls vs. multiple sclerosis], and (iii) the interaction effect of time by disease. The interaction term was added to investigate the temporal trajectories of the dependent measures in the two study groups and to test the hypothesis that multiple sclerosis patients have an accelerated increase in magnetic susceptibility compared to controls. In addition, specific within-study-group time effects were tested using the paired t-test, and between-study-group baseline and follow-up differences were tested using the Student's t-test. To test whether iron accumulation occurs predominantly in the earlier disease stages (and subsequently slows down) (Hagemeier et al., 2012b, Al-Radaideh et al., 2013), we repeated all susceptibility analyses separately for RR- and SPMS sub-groups. Effect sizes were estimated using Hedges' g. Associations between cross-sectional and longitudinal MRI and clinical measures were explored using Pearson or Spearman correlation coefficients, depending on data distribution. P-values were corrected for multiple comparisons using the false discovery rate method (then denoted as q-values) (Benjamini and Hochberg, 1995) and results were considered statistically significant when q < 0.05; uncorrected p < 0.05 was considered a trend.

Results

Clinical and demographic measures at baseline

Table 1 shows baseline clinical, demographic, and MRI differences between controls and multiple sclerosis patients. The study groups were similar in age and distribution of sex (p > 0.1), with a mean age of 43.7 ± 12.3 years for controls and 44.2 ± 10.2 years for patients. Twenty-four (24) controls (60%) and 81 (67.5%) multiple sclerosis patients were female. Mean disease duration of the multiple sclerosis group was 12.8 ± 9.4 years. Mean follow-up time was 1.9 ± 1.2 years for controls and 2.1 ± 1.1 years for multiple sclerosis patients (p > 0.1).

Table 1

Clinical, demographic, and volumetric MRI data at baseline and change over follow-up.

	Controls	Multiple sclerosis	p-value
N	40	120	NA
Female, n (%)	24 (60%)	81 (67.5%)	0.387c
Age	43.7 ± 12.3	44.2 ± 10.2	0.819d
RRMS, n (%)	NA	98 (81.7%)b	NA
Baseline EDSS, median ± IQR	NA	2.5 ± (1.5–4.0)	NA
Follow-up EDSS, median ± IQR	NA	3.0 ± (1.5–4.0)	NA
Disease duration, years	NA	12.8 ± 9.4	NA
Baseline relapses in the last year, median ± IQR	NA	0 ± (0–1)	NA
Follow-up relapses in the last year, median ± IQR	NA	0 ± (0–0)	NA
Treated with DMT, n (%)a	NA	103 (85.3%)	NA
Follow-up time, years	1.9 ± 1.2	2.1 ± 1.08	0.516d
Baseline T2 lesion volume, ml	0.3 ± 0.7	13.9 ± 18.1	< 0.001d, e
Follow-up T2 lesion volume, ml	0.7 ± 1.7	12.6 ± 14.8	< 0.001d, e
Whole brain volume, ml	1532.0 ± 76.2	1457.8 ± 91.4	< 0.001d, e
Gray matter volume, ml	768.0 ± 49.2	743.6 ± 54.6	0.048d, e
White matter volume, ml	770.1 ± 43.8	723.7 ± 52.5	< 0.001d, e
Lateral ventricle volume, ml	33.6 ± 16.2	48.1 ± 20.9	< 0.001d, e
PB-VC	− 0.71 (0.82)	− 1.12 (1.16)	0.139d, e
PGM-VC	− 1.29 (2.04)	− 1.45 (1.66)	0.807d, e
PWM-VC	0.29 (1.54)	0.50 (2.07)	0.764d, e
PV-VC	3.74 (6.11)	5.12 (7.40)	0.535d, e

Abbreviations: NA = not applicable; RRMS = relapsing-remitting multiple sclerosis; EDSS = Expanded Disability Status Scale; IQR = inter-quartile range; DMT = disease modifying therapy; PB-VC = percent brain volume change, PV-VC = percent ventricle volume change; PGM-VC = percent gray matter volume change; PWM-VC = percent white matter volume change.

Results are based on baseline and percent change data, and are presented as mean ± SD unless otherwise noted. Brain volumes were normalized for head size.

Patients were taking the following disease modifying treatments: interferon beta: 40, natalizumab: 21, glatiramer acetate: 34, combination: 3, other: 5, no therapy: 17. 11 (9.2%) had a change in DMT from baseline to follow-up.

Patients not identified as RRMS were secondary progressive (SP).

Chi-squared test.

Independent Student t-test.

MRI p-values were corrected using false discovery rate resulting in q-values.

As expected, SP patients were older than RR patients (52.7 ± 6.7 vs. 42.0 ± 9.6 years; p < 0.001). Although non-significant (p = 0.21), there were more female (70%, n = 69) RRMS than SPMS patients (54.5%, n = 12).

Cross-sectional analysis

Magnetic susceptibility in DGM nuclei of controls correlated highly with published histochemical brain iron concentrations (R2 > 0.87, p < 0.001, Supplementary Fig. 1) (Hallgren and Sourander, 1958), illustrating that DGM susceptibility is dominated by iron in the normal brain and that QSM is a sensitive approach to study tissue iron.

Supplementary Fig. 1

Magnetic susceptibility of healthy controls plotted against estimated iron concentrations based on analytical equations given in (Hallgren and Sourander, 1958) derived from histochemical iron quantification. Thalamus iron concentration was estimated from Figure 8 in (Hallgren and Sourander, 1958), as described in Supplement 2e. Blue squares: thalamus; green triangles: caudate; tan circles: putamen; and purple crosses: globus pallidus. Dotted line represents the best fit (R2 = 0.91: y = − 26.1 ppb + 6.2 ppb ∙ 100 g/mg ∙ [Fe], where [Fe] is the iron concentration in mg/100 g).

Table 2 presents an overview of cross-sectional magnetic susceptibilities for controls and multiple sclerosis patients. At both time points, susceptibilities of the caudate (g ≥ 0.54, q ≤ 0.01) and globus pallidus (g ≥ 0.92, q < 0.001) were significantly higher, and susceptibility of the thalamus was significantly lower (g ≥ 0.71, q < 0.001) in multiple sclerosis compared to controls. Similar results were obtained when restricting the multiple sclerosis group to RRMS patients vs. the more similarly aged controls (data not shown).

Table 2

Baseline and follow-up magnetic susceptibility in controls and multiple sclerosis.

		Baseline	Follow-up	Difference (effect size)	qa	qc	qd
Thalamus	Controls	− 0.4 ± 6.7	0.2 ± 6.6	+ 0.6 (g = 0.09)	0.603	0.877	0.297
	Multiple sclerosis	− 6.1 ± 8.3	− 6.9 ± 7.8	− 0.8 (g = 0.10)	0.251
	Difference (effect size)	− 5.7 (g = 0.71)	− 7.1 (g = 0.94)
	qb	< 0.001	< 0.001
	qe	< 0.001
Caudate	Controls	36.4 ± 8.8	38.9 ± 9.9	+ 2.5 (g = 0.27)	0.045	< 0.001	0.628
	Multiple sclerosis	44.3 ± 12.8	45.9 ± 13.4	+ 1.6 (g = 0.12)	0.014⁎
	Difference (effect size)	+ 7.9 (g = 0.66)	+ 6.9 (g = 0.54)
	qb	< 0.001	0.014⁎
	qe	0.001
Putamen	Controls	48.7 ± 11.8	49.2 ± 11.9	+ 0.5 (g = 0.04)	0.737	0.347	0.807
	Multiple sclerosis	51.8 ± 15.6	52.8 ± 16.2	+ 1.0 (g = 0.06)	0.224
	Difference (effect size)	+ 3.1 (g = 0.21)	+ 3.6 (g = 0.23)
	qb	0.469	0.403
	qe	0.424
Globus pallidus	Controls	103.3 ± 15.7	102.6 ± 13.6	− 0.7 (g = 0.05)	0.807	0.877	0.611
	Multiple sclerosis	124.2 ± 24.5	125.4 ± 25.5	+ 1.2 (g = 0.05)	0.541
	Difference (effect size)	+ 20.9 (g = 0.92)	+ 22.8 (g = 0.98)
	qb	< 0.001	< 0.001
	qe	< 0.001

Note: P-values were corrected using false discovery rate resulting in q-values. Between-subject q-values are presented below descriptive results, while within-subject q-values are presented to the right of descriptive results. Results described as mean ± SD, susceptibility is stated in ppb. Effect sizes were computed using Hedge's g.

q-value of the within group difference between baseline and follow-up (paired t-test).

q-value of the difference between controls and multiple sclerosis (independent t-test).

Main effect of time (mixed factorial ANOVA).

q-value of the interaction effect of time by disease (mixed factorial ANOVA).

Main effect of disease (mixed factorial ANOVA).

Statistical significance of caudate magnetic susceptibility among multiple sclerosis patients was lost (p > 0.05) when adjusting for caudate volume. Other structures were unaffected.

Volumes of all DGM regions were reduced in multiple sclerosis compared to controls at both baseline and follow-up (between − 5.3% and − 10.5%, g ≥ 0.41, q ≤ 0.05). Detailed results of the cross-sectional analysis of DGM and lesion volumes may be found in Supplement 1.

Intra-subject temporal evolution

Intra-subject 2-year changes of susceptibility (Table 2 and Fig. 2) reached significance only in the caudate of controls (+ 2.5 ppb, g = 0.27, q < 0.05) and multiple sclerosis (+ 1.6 ppb, g = 0.12, q = 0.01). The significance of temporal changes of caudal susceptibility was lost within multiple sclerosis patients when the model was adjusted for caudate volume loss, suggesting an association between atrophy and susceptibility increase in this structure (Table 2). No statistically significant interaction effects of susceptibility were observed in any region, indicating a similar rate of increase in both groups over time. Two-year changes of DGM volumes were significant among MS but not HC. Details of the temporal evolution analysis of DGM- and lesion- volumes are presented in Fig. 3 and Supplement 1.

Fig. 2

Susceptibility change trajectories as a function of age for the thalamus, caudate, putamen, and globus pallidus. The left column depicts controls while the right column depicts multiple sclerosis patients. Gray circles indicate baseline; black circles indicate follow-up. The dashed linear fits represent the mean change in susceptibility as a function of age. Regression coefficients are as follows: controls: thalamus = − 0.107 ppb/years ∙ [age] + 4.713 ppb, caudate = 0.363 ppb/years ∙ [age] + 21.305 ppb, putamen = 0.501 ppb/years ∙ [age] + 26.310 ppb, globus pallidus = 0.416 ppb/years ∙ [age] + 84.117 ppb; multiple sclerosis: thalamus = − 0.227 ppb/years ∙ [age] + 3.834 ppb, caudate = 0.480 ppb/years ∙ [age] + 23.157 ppb, putamen = 0.501 ppb/years ∙ [age] + 29.447 ppb, globus pallidus = 0.629 ppb/years ∙ [age] + 96.024 ppb.

Fig. 3

Trajectories of individual normalized volumes as a function of age for the thalamus, caudate, putamen, and globus pallidus. The left column depicts controls while the right-hand side depicts multiple sclerosis patients. Gray circles indicate baseline; black circles indicate follow-up. The dashed linear fits represent the mean change in volume as a function of age. Regression coefficients are as follows: controls: thalamus = − 0.066 ml/years ∙ [age] + 23.647 ml, caudate = − 0.037 ml/years ∙ [age] + 10.738 ml, putamen = − 0.036 ml/years ∙ [age] + 14.46 ml, globus pallidus = 0.005 ml/years ∙ [age] + 4.425 ml; multiple sclerosis: thalamus = − 0.080 ml/years ∙ [age] + 22.270 ml, caudate = − 0.041 ml/years ∙ [age] + 10.266 ml, putamen = − 0.044 ml/years ∙ [age] + 13.991 ml, globus pallidus = − 0.001 ml/years ∙ [age] + 4.291 ml.

Similar to the total multiple sclerosis group, the caudate showed significantly increased susceptibility over time in RRMS (baseline: 43.3 ± 13.5 ppb, follow-up: 44.9 ± 13.9 ppb, + 1.6 ppb, g = 0.12, q = 0.049) and a tendency (smaller group) for increased susceptibility in SPMS patients (baseline: 48.7 ± 8.4 ppb, follow-up: 50.6 ± 9.7 ppm, + 1.9 ppb, g = 0.21, p = 0.092; Supplement Table 2). However, different from the total patients group, in RRMS the putamen also showed a tendency toward higher susceptibility at follow-up compared to baseline (baseline: 50.5 ± 15.6 ppb, follow-up: 52.0 ± 16.3 ppb, + 1.5 ppb, g = 0.09, p = 0.017).

Associations of susceptibility with clinical and other MRI outcomes

In multiple sclerosis, changes of EDSS over two years did not correlate with changes of susceptibility in any of the regions. No statistically significant correlations were found in multiple sclerosis between susceptibility and disease duration nor percent ventricle volume change, Gd lesion number change, and any of GM, WM, neocortical, brain parenchymal, lateral ventricle, putamen, globus pallidus, T1, and T2 lesion volumes. A significant correlation was found between susceptibility and volume change in the caudate (r = − 0.247, q = 0.04). There was a trend for SIENA-determined percent brain volume loss with decreased thalamic susceptibility (r = 0.208, p = 0.027, q = 0.09). Among controls, within-structure change of volume and susceptibility over two years was correlated in the putamen (r = − 0.477, q < 0.01) as well as caudate (r = − 0.247, p < 0.01). Within-structure associations did not reach significance in other DGM structures.

Discussion

The present work reports the first longitudinal measurement of magnetic susceptibility in the human brain. Contrary to most other imaging-based quantities, magnetic susceptibility is a physical property that is invariant to experimental parameters such as sequence parameters or field strength. Although QSM methodologies differ between labs, results obtained in different experimental settings may be compared with one another (Deh et al., 2015, Lin et al., 2015, Santin et al., 2016). Furthermore, the addition of QSM to the quantitative neuroimaging toolbox sheds light on the underlying cause of previously reported findings with R2*, because of the complementary dependence of QSM and R2* on myelin and calcium. While iron increases both measures, myelin and calcium decrease susceptibility and increase R2*. If both iron and myelin/calcium concentrations change, the relative amount of each substance determines the direction of changes of the total voxel values.

Detection of iron concentration changes in the normal brain over two years

Temporal susceptibility changes reached statistical significance only in the caudate (+ 2.5 ppb, 95% CI [0.55, 4.45] ppb). Hallgren and Sourander's landmark brain iron study (Hallgren and Sourander, 1958) allows for an estimation of the expected susceptibility increase in deep gray matter nuclei due to normal aging between 44 and 46 years of age (the mean age of our cohort at the two time points). For the caudate, the Hallgren and Sourander study predicts + 0.71 ppb (Supplement 2a). The average susceptibility difference observed in the present study is not inconsistent with this value (in terms of 95% statistical confidence). Moreover, a recent meta-analysis of 20 MRI studies of age-related iron (Daugherty and Raz, 2013) found that in vivo estimated differences in the caudate were generally larger than estimates reported by Hallgren and Sourander (Daugherty and Raz, 2015), indicating that their measures underestimated the iron concentration.

An increase of the iron concentration may generally be caused by an actual deposition of iron (e.g. more ferritin or higher ferritin loading factor) or by a removal of structures that do not contain iron, resulting in an increased density of the existing iron. Considering that we did not detect caudal atrophy over two years in controls, it may be concluded that the increased iron concentration in this region is in fact caused by iron deposition. Daugherty et al. (Daugherty and Raz, 2016) did not find longitudinal changes in the caudate, but did in the putamen over seven years in an older cohort (70.85 ± 9.91 years at end of study) using R2*. However, they also found significant atrophy in both regions, indicating loss of myelin, which has a decreasing effect on R2*, rendering the detection of simultaneous iron accumulation and myelin loss difficult.

Temporal susceptibility changes between baseline and follow-up did not reach significance in other regions as our study was not sufficiently powered to detect normal iron accumulations (see Supplement 2b), and the two-year time frame may not be sufficient to reliably detect susceptibility changes due to changes being too small relative to previously reported QSM reproducibility (Deh et al., 2015, Lin et al., 2015).

Increased iron concentration in caudate and globus pallidus of multiple sclerosis patients

In line with previous cross-sectional QSM studies (Al-Radaideh et al., 2013, Langkammer et al., 2013, Rudko et al., 2014), we found increased caudal and pallidal magnetic susceptibilities in multiple sclerosis patients compared to controls at both time-points. The susceptibility difference of the caudate (at baseline) was smaller (+ 7.9 ppb, 95% CI [3.60, 12.20] ppb) than that reported by Langkammer et al. (+ 15.0 ppb, 95% CI [8.05, 21.98] ppb). The difference in pallidal susceptibility was higher in our study compared to Rudko et al. (+ 20.9 ppb, 95% CI [12.73, 29.07] ppb vs. + 13.0 ppb, 95% CI [10.5, 15.5] ppb). Overlapping caudal and pallidal confidence intervals between these studies and our work suggest that results are consistent. Although, direct comparisons with results by Rudko et al. (2014), Cobzas et al. (2015), and Al-Radaideh et al. (2013) are challenging, as these authors referenced susceptibility to white matter and internal capsule (regions known to be affected by iron depletion and myelin loss in multiple sclerosis and by head orientation effects due to susceptibility anisotropy in WM).

Increased magnetic susceptibility can have direct or indirect biophysical causes: due to the diamagnetism of myelin, a reduction of myelin and/or a deposition of (paramagnetic) iron directly increase susceptibility. Marked demyelination in the caudate nucleus of multiple sclerosis patients has recently been reported post mortem (Haider et al., 2014), although the generally low myelin content of the DGM nuclei suggests that susceptibility changes in these regions are predominantly related to iron. Also physical shrinkage of a structure (unrelated to change in myelin and not affecting iron-containing cells) could theoretically increase magnetic susceptibility indirectly by increasing the concentration of existing tissue iron (no actual iron deposition). An increased iron concentration has been speculated to amplify iron toxicity, potentially causing a vicious cycle of further structural shrinkage.

Increasing susceptibility of caudate is likely a primary consequence of atrophy and not due to iron deposition

We observed the highest atrophy rate over two years (≥ 1% per year) in the caudate of patients, followed by the globus pallidus (regions that also showed significant cross-sectional susceptibility differences) and putamen. Assuming that the total amount of iron remained constant, this physical shrinkage alone can explain average susceptibility increases of + 1.5 ppb, + 3.7 ppb, and + 1.2 ppb in the caudate, globus pallidus, and putamen, respectively, in patients (Supplement 2c). Although only the caudate reached statistical significance, the observed susceptibility increases over two years (+ 1.6 ppb, + 1.0 ppb, and + 1.2 ppb, respectively) were similar to those estimations, suggesting that measured susceptibility increases in these regions are dominated by physical condensation of iron. However, histopathology showed that the caudate suffers from substantial demyelination with reduced iron density in DGM lesions (Haider et al., 2014), suggesting that the observed susceptibility increase is composed of contributions from both increased density of existing iron (due to physical shrinkage) and partially compensating components of myelin loss and iron depletion. In particular, the idea of iron depletion is supported by our finding of a significantly higher rate of atrophy in patients, coupled with a similar rate of increase of susceptibility compared to controls. Considering that iron is also deposited in the caudate during normal aging (see above), atrophy and myelin loss alone would result in a considerably higher susceptibility increase compared to controls. Further support is provided by the absence of caudal changes in previous R2*-based longitudinal studies (e.g. (Walsh et al., 2014, Khalil et al., 2015, Daugherty and Raz, 2016))—the opposing effect of myelin loss and iron on relaxometry renders the detection of such processes difficult with R2*. Overall, our hypothesis that multiple sclerosis patients have an accelerated increase in magnetic susceptibility compared to controls could not be confirmed.

Caudate is likely a continuous target throughout the disease course

In line with previous studies (Khalil et al., 2009, Burgetova et al., 2010, Hagemeier et al., 2012b, Hagemeier et al., 2013c), we found significant temporal volume changes in all DGM structures, indicating that pathological processes are ongoing in our sample of multiple sclerosis patients. Moreover, the statistically significant within-structure correlation between susceptibility and atrophy in the caudate further corroborates this association and suggests that caudate pathology is a continuous target throughout the multiple sclerosis disease course. Interestingly, an estimation of the caudate susceptibility increase over two years from the baseline cross-sectional susceptibility difference between study groups (+ 1.2 ppb; Supplement 2d) yields a value that is similar to the measured caudate susceptibility increase over two years (+ 1.5 ppb). However, considering that a similar rate of increase of susceptibility was observed in the caudate of patients and controls, it may be concluded that while pathology is ongoing, the observed caudal cross-sectional susceptibility difference must have been primarily caused at earlier stages of the disease.

Disease-related susceptibility changes in globus pallidus and putamen may be masked by an age-dependent decrease of the atrophy rate in patients

Although correlations and temporal changes did not reach statistical significance outside the caudate, it should not be concluded that the observed (ongoing) atrophy in these regions is not associated with increased iron concentrations. For example, cross-sectionally higher susceptibility in the putamen has recently been reported in younger (and smaller) multiple sclerosis cohorts (Langkammer et al., 2013, Cobzas et al., 2015). In our study, putamen did not reach significance, but there was a tendency of a difference between relapsing remitting and secondary progressive patients (Supplement Table 2). Furthermore, significantly increased longitudinal R2* has recently been found (over two years) in the globus pallidus of clinically isolated syndrome patients (Khalil et al., 2015) and multiple sclerosis patients with a mean disease duration of 5.8 years (Walsh et al., 2014), but not in patients with (a higher) mean disease duration of 7.5 years (over three years) (Khalil et al., 2015). These observations suggest that disease duration or age-dependent effects mask temporal susceptibility changes at later stages of the disease.

The effect sizes of both temporal susceptibility and volume changes were substantially lower in the putamen (g = 0.06 and g = 0.16, respectively) and globus pallidus (g = 0.05 and g = 0.18, respectively) compared to the caudate (g = 0.12 and g = 0.23). Increasing variation of iron concentrations in the DGM with age, which could explain the small effect sizes, is known from Hallgren and Sourander's study in healthy subjects (Hallgren and Sourander, 1958) and has been confirmed with R2* (Daugherty et al., 2015) and QSM (Li et al., 2014). However, while it is clear that atrophy-related changes in the thalamus of patients remained relatively low (Supplement Table 1), a potentially reduced disease-related atrophy rate at later ages may also be responsible for undetected susceptibility changes in the globus pallidus and putamen. This hypothesis of a slowdown of iron accumulation (potentially due to atrophy-related condensation of existing iron) at later stages of the disease is supported by a recent neuropathological study reporting that iron loading increased with age among controls, but not multiple sclerosis patients (Haider et al., 2014). Future longitudinal studies contrasting different disease stages and deep gray matter structures as well as a meta-analysis of the current literature on brain iron in multiple sclerosis should shed light on this conundrum.

In the globus pallidus and putamen, where atrophy rates were relatively high, susceptibility changes may have remained undetected due to moderately high inter-subject variations of susceptibility values (globus pallidus; Table 2) and volumes (putamen; Supplement Table 1). Although, standard deviations of susceptibility values in these two regions (approximately 16 ppb and 25 ppb) were similar to those reported in younger groups (Langkammer et al., 2013, Rudko et al., 2014), additionally, the range of values was consistent with a previous report in controls (Li et al., 2014).

When separately analyzing RRMS and SPMS (and thereby splitting our cohort into more homogeneous younger and older groups) the putamen showed increasing susceptibility over time in RRMS (but not SPMS), while the caudate showed an increasing susceptibility in both RRMS and SPMS. It has been suggested that iron increases may already be pronounced in early disease stages (Hagemeier et al., 2012b, Al-Radaideh et al., 2013). In this context, the current observations suggest that temporal changes are similar between RRMS and SPMS in the caudate, more pronounced in the putamen among RRMS (earlier disease stage), and are stable in the thalamus and pallidus. This confirms that the caudate is a continuous target in disease, while putamen susceptibility change may drop to a rate similar to the expected rate in normal aging. Furthermore, the older SPMS group did not present with higher standard deviations in our sample (data not shown), suggesting that although the sample may be smaller, the SPMS measurements are representative.

Decreased thalamic magnetic susceptibility in multiple sclerosis

The literature is inconsistent regarding findings of thalamic tissue changes. Both Al-Radaideh et al. (2013) and Langkammer et al. (2013) did not find significantly different (global) thalamus susceptibility in clinically isolated syndrome and multiple sclerosis patients, respectively, whereas Rudko et al. (2014) reported substantially increased thalamic susceptibility (+ 11 ppb) in a group of multiple sclerosis patients. Furthermore, significantly increased thalamic R2* has been found in patients with multiple sclerosis (Langkammer et al., 2013) but decreasing thalamic R2* values over three years have been reported in a mixed group of patients with clinically isolated syndrome and multiple sclerosis (Khalil et al., 2015). A coherent explanation of these seemingly contradictory findings has not yet been attempted and is difficult because iron concentration in the thalamus is known to follow a peculiar non-linear aging trajectory in normal controls (Hallgren and Sourander, 1958, Hagemeier et al., 2013a). In addition, published studies relied on groups with different demographic and disease characteristics. In particular, the groups studied in the literature were substantially younger than our group (average age between 30 and 40 years, compared to 44 years).

Disregarding group differences and using the traditional myelin-iron model to explain the changes of susceptibility and R2* reported in the literature lead to the conclusion of simultaneous iron accumulation and myelin loss in the thalamus. This is because increased susceptibility is indicative for iron accumulation and/or myelin loss, increased R2* is indicative for iron accumulation, and decreasing R2* is indicative for myelin loss. In this traditional model, our finding of lower thalamic susceptibility in patients would indicate a lower iron concentration and/or a higher myelin density, but not myelin loss (which would increase the susceptibility).

Neuronal loss in the thalamus of multiple sclerosis patients has been described in neuropathological studies (Cifelli et al., 2002, Vercellino et al., 2009). Gray matter pathology should be associated with an increasing ratio of white to gray matter volume (increased myelin density) and may, in addition, involve direct removal of iron-containing cells; both effects would explain the reduced magnetic susceptibility in our multiple sclerosis group. A lower iron concentration would result in decreased R2*, which is consistent with the report of decreasing R2* in patients by Khalil et al. (2015). Iron depletion combined with increased myelin density due to gray matter loss may be the main culprit of reduced thalamic susceptibility.

Another at this point arguably speculative explanation of decreased susceptibility, is via the deposition of insoluble calcium phosphate—similar to myelin, calcium phosphate decreases susceptibility (diamagnetic) (Schweser et al., 2010) and increases R2*. Excessive axonal influx of Ca2 + has been associated with a cascade of compensatory mechanisms due to reduced signal transduction following demyelination in multiple sclerosis (Trapp and Stys, 2009). Mitochondria seem to play a critical role in the conversion of free calcium cations into insoluble matrix calcium phosphate, because increased cytosolic calcium concentrations can trigger a mitochondrial permeability transition, which has been shown to result in the formation of calcium phosphate (Panov et al., 2004). Supporting evidence for this hypothesis is given by the histopathological finding of focal demyelinated lesions in the thalamus of multiple sclerosis patients (Haider et al., 2014) and secondary calcium influx into the thalamus observed in a model of traumatic brain injury even if the primary injury occurred elsewhere (Osteen et al., 2001), suggesting that the thalamus is particularly susceptible for calcium influx following neuronal injury.

To further investigate whether iron depletion and/or calcium deposition can explain the reported thalamic susceptibility and R2* findings in the different studies, additional analyses were carried out that are reported in the Supplementary material.

Limitations of the study

The comparison of our results with the literature (Supplementary Fig. 2) shows that it is important to study iron in homogeneous groups of patients with multiple sclerosis at different disease stages, which has not been done to date mostly due to the challenge of collecting sufficiently large samples. In particular, our Monte Carlo simulation (Supplementary material) revealed how much group differences (with a wide spread of ages) can differ from individual trajectories (Supplementary Fig. 3), motivating studies with narrower age ranges and longer follow-up times. Future studies with longer follow-up periods could also increase power and decrease spread in order to better detect longitudinal changes between groups. This need is apparent due to longitudinal changes in the present study being subtle in nature (within or below the range of inter and intra-scanner QSM reproducibility (Deh et al., 2015, Lin et al., 2015, Santin et al., 2016).

Supplementary Fig. 2

Monte-Carlo simulation of expected age-dependent group-average trajectories of susceptibility (a, c) and R2* (b, d) in the thalamus. The iron depletion model is shown in the top row (a, b), the calcium model in the bottom row (c, d). Straight red lines represent the simulated cross-sectional group-control differences. Straight and dashed black lines are the simulated normal and MS group values. Data points with error bars (95% CI; estimated from reported standard deviations/errors) represent values found at both time points in the present study and previously reported cross-sectional differences in Al-Radaideh et al. (2013) (A), Langkammer et al. (2013) (B), Rudko et al. (2014) (C), and Khalil et al. (2015) (D). Details of the simulations are reported in Supplement 3.

Another limitation is the ROI-based analysis employed in this and most other studies. The employed FSL FIRST method has been shown to suffer from inaccuracies in the presence of severe brain abnormalities, such as atrophy (Amann et al., 2015), which may introduce disease-related anatomical bias. Inaccuracies of the method stem primarily from a poor transformation to NMI space (Feng et al., 2017), which suggests that the inaccuracies manifest randomly, not systematically. Also, mean image intensity (susceptibility) measures within ROIs, are likely to be substantially less affected than volumetric measures. While these deviations may have masked group differences in our study, most previous studies employed the same methodology and reported similar standard deviations.

A general limitation of the a priori definition of anatomical regions is that ROIs may not reflect the actual tissue pathology, which may be restricted to certain anatomical sub-regions. For example, Al-Radaideh et al. (2013) did not find global thalamus differences, but reported higher susceptibility in the pulvinar of the thalamus (+ 13 ppb, d = 0.86) in clinically isolated syndrome patients. However, automated segmentation of sub-regions with individual DGM structures remains challenging. Future voxel-based analyses may allow us to study the spatial distribution of changes in more detail.

In addition, while the thalamic results are intriguing, we cannot rule out that results were influenced by thalamic and DGM lesions, which are known to occur in multiple sclerosis (Haider et al., 2014). That said, we did not observe thalamic lesions on QSM.

In this study, susceptibility values were referenced to the whole brain average susceptibility. The impact of referencing on the study results is subject to an active debate in the field. Some authors have recommended the use of cerebrospinal fluid (CSF) or the internal capsule (Straub et al., 2016), while others recommended frontal deep white matter (Deistung et al., 2013). The choice of reference region has substantial impact on the resulting group differences because pathology-related changes in the reference region propagate into differences in referenced susceptibility values. For example, an increased reference susceptibility due to pathology-related demyelination would decrease the apparent susceptibility of referenced regions in the patient group, potentially masking iron-related susceptibility increase. In addition, inter-subject variability of the absolute reference susceptibility value increases the inter-subject variability of referenced values. Al-Radaideh et al. (2013) and Cobzas et al. (2015) used the internal capsule as a reference region. However, apart from potential demyelination in this white matter region and the proximity of the internal capsule to the globus pallidus (increased susceptibility in multiple sclerosis) requiring accurate segmentation, the susceptibility in this highly anisotropic region may depend on head orientation in the magnet. Furthermore, due to the limited size of the region the mean value calculation is prone to localized reconstruction artifacts. We presumed referencing to the whole brain average represents the best trade-off between susceptibility anisotropy effects, localized artifacts, and pathology related susceptibility changes. However, global white matter demyelination could still influence the region.

The reasons for change of iron concentration in multiple sclerosis are still not well understood. Our results indicate that the disturbance of iron homeostasis in multiple sclerosis is a complex phenomenon and findings with iron-sensitive MRI may not necessarily reflect only deposition of iron. Iron concentrations may be changing mostly during pathological events such as in the chronic pro-inflammatory phase of white matter lesions (Mehta et al., 2013) or in the early stages of structural atrophy (Hagemeier et al., 2012b, Hagemeier et al., 2013c, Blazejewska et al., 2015), perhaps both as an instigator of oxidative damage as well as a byproduct of (myelin) destruction. Although further research will have to investigate the exact mechanisms at play in multiple sclerosis, the results of the present study corroborate that structural atrophy is intricately linked to increased iron density.

Conclusions

We hypothesized a presence of higher susceptibility in multiple sclerosis patients and that the susceptibility gap was augmented over time. Our results partially support this notion: while susceptibility was found to be cross-sectionally higher in the caudate and globus pallidus, consistent with previous reports, we found lower susceptibility in the thalamus. Results indicate continued changes in the caudate. The increased susceptibility in the caudate may be attributed to a physical increase of iron density due to structural atrophy rather than actual deposition of iron. The comparison of atrophy and susceptibility changes in patients supports previous findings of additional myelin loss in this region. Thalamic findings indicate a multi-phase pathology that can be explained by simultaneous myelin loss, iron accumulation and either iron depletion (e.g. due to gray matter atrophy) or calcium deposition (e.g. as a secondary effect of demyelination) at later stages of the disease. Further clarification necessitates future investigation.

The following are the supplementary data related to this article.

Supplementary material

Detailed results of the Monte Carlo simulation of susceptibility and R2* changes in the thalamus for the iron depletion (top row) and the calcium deposition (bottom row) hypotheses, respectively. (a) Visualization of the events and their relative strength with iron in black, demyelination in blue, and calcium phosphate in red. Shown is the change (derivative) of the tissue composition changes (e.g. iron concentration follows sigmoid function), reflecting the strength of the pathological event. (b) Individual contributions of iron, myelin, and calcium to the average voxel susceptibility (see panel c) in normal subjects (straight lines) and multiple sclerosis patients (dashed lines). (c) Average voxel susceptibility predicted by the models for normal subjects (straight black) and patients (dashed black), and difference due to the disease (straight red). Panels (d) and (e) are analogous plots to (b) and (c) for R2*. The trajectories in panels (c) and (e; difference shifted by + 17/s for improved visualization) were the basis for the Monte Carlo simulation, which is shown in Supplementary Fig. 2.