Diffusion Kurtosis Imaging maps neural damage in the EAE model of multiple sclerosis.

Diffusion kurtosis imaging (DKI) is an imaging modality that yields novel disease biomarkers and in combination with nervous tissue modeling, provides access to microstructural parameters. Recently, DKI and subsequent estimation of microstructural model parameters has been used for assessment of tissue changes in neurodegenerative diseases and associated animal models. In this study, mouse spinal cords from the experimental autoimmune encephalomyelitis (EAE) model of multiple sclerosis (MS) were investigated for the first time using DKI in combination with biophysical modeling to study the relationship between microstructural metrics and degree of animal dysfunction. Thirteen spinal cords were extracted from animals with varied grades of disability and scanned in a high-field MRI scanner along with five control specimen. Diffusion weighted data were acquired together with high resolution T2* images. Diffusion data were fit to estimate diffusion and kurtosis tensors and white matter modeling parameters, which were all used for subsequent statistical analysis using a linear mixed effects model. T2* images were used to delineate focal demyelination/inflammation. Our results reveal a strong relationship between disability and measured microstructural parameters in normal appearing white matter and gray matter. Relationships between disability and mean of the kurtosis tensor, radial kurtosis, radial diffusivity were similar to what has been found in other hypomyelinating MS models, and in patients. However, the changes in biophysical modeling parameters and in particular in extra-axonal axial diffusivity were clearly different from previous studies employing other animal models of MS. In conclusion, our data suggest that DKI and microstructural modeling can provide a unique contrast capable of detecting EAE-specific changes correlating with clinical disability.

Introduction

Multiple sclerosis (MS) is a demyelinating, inflammatory, neurodegenerative disease of the human central nervous system (CNS) affecting millions of people worldwide. The pathophysiology of MS is often complex, and involves, among other factors, myelin loss, axonal damage, appearance of transient or permanent lesions, and brain atrophy. Effective treatment of MS is still lacking (Compston and Coles, 2002), although a range of disease-modifying therapies have been introduced (Berger, 2011; Noyes and Weinstock-Guttman, 2013). These therapies are based on immunomodulatory, anti-inflammatory, and immunosuppressive drugs. The success of such treatments depends on early diagnosis and careful monitoring of the patient.

A range of MS animal models characterized by different mechanisms of induction and pathology (Lassmann and Bradl, 2016) have been developed to overcome the limitations of clinical tissue assessment. Experimental autoimmune encephalomyelitis (EAE) is a group of the most compelling and commonly used animal MS models (Baker and Amor, 2014; Kipp et al., 2016; Lassmann and Bradl, 2016). Unlike other animal models, in addition to inflammatory lesions and demyelination, EAE includes axonal damage (Bergers et al., 2002; Kipp et al., 2016) which is one of the hallmarks of MS. Therefore, using EAE to assess MS biomarkers can provide unique insights into the MS pathology.

Due to its noninvasiveness and ability to contrast soft tissues, Magnetic Resonance Imaging (MRI) is extensively used for diagnosis and monitoring of MS (Bakshi et al., 2008; Polman et al., 2011). Standard T1- or T2-weighted MRI images are capable of revealing brain atrophy and lesions, which are heterogeneous areas harboring demyelination, inflammation, gliosis and axonal injury (Filippi et al., 2012; Inglese and Bester, 2010). However, diffuse microstructural changes outside the T1 or T2 intensity lesions in gray matter (GM) and so-called normal appearing white matter (NAWM) (Allen et al., 2001) have been observed with histology. Studies (De Stefano et al., 2006; Kipp et al., 2016; Miller et al., 2003) showed that the diffuse damage in NAWM and GM contributes to disability accumulation and chronic disease progression.

Diffusion weighted imaging (DWI) can provide quantitative microstructural information by sensitizing MRI signals to the displacement of water molecules. The underlying signal attenuation is often approximated by a Gaussian distribution, which forms the basis of diffusion tensor imaging (DTI) (Basser and Pierpaoli, 1996). While DTI metrics are widely used, the mentioned Gaussian approximation is valid only in a limited regime of low diffusion weighting. At higher gradient strengths (or b-values), tissue microstructure and compartmentalization cause deviations from Gaussianity. These deviations are utilized by a framework known as diffusion kurtosis imaging (DKI) (Jensen et al., 2005; Jensen and Helpern, 2010). Combined with tissue modeling, DKI provides access to microstructural parameters and yields novel disease biomarkers.

DKI based biomarkers have been shown to improve the diagnostic assessment in a range of neurological disorders (Delgado y Palacios et al., 2014; Grossman et al., 2011; Khan, 2016; Surova et al., 2016; Tietze et al., 2015; Wang et al., 2011). In MS, they improved characterization of GM and NAWM damage (Raz et al., 2013; Yoshida et al., 2013) and correlated with cognitive impairment (Bester et al., 2015). In MS animal models, DKI biomarkers were associated with chronic injury (Falangola et al., 2014; Guglielmetti et al., 2016; Jelescu et al., 2016) and neurite myelin content (Kelm et al., 2016). However, even though (Wu and Cheung, 2010) employed EAE to show that DKI is able to enhance lesion detection and other DWI methods revealed pathological changes in EAE (Biton et al., 2005; Budde et al., 2009), DKI and white matter (WM) models have never been used to investigate EAE-induced disability.

In this study, we hypothesized that novel metrics obtained using DKI could predict the EAE disability. If this hypothesis is found to be true, and provided our methods are successfully translated to clinic setting, the results of this study can support usage of DKI-derived biomarkers as a diagnostic tool for MS. A detailed description of the chosen DKI-derived metrics is provided in the next chapter following the definition of the metrics and a description of their relationship to the pathology.

The ability of DKI to provide quantitative biomarkers of dysfunction in EAE model of MS was investigated by testing the correlation between the biomarkers and behavioral markers of disease severity. Inside lesions no biomarkers showed correlation to disability. In NAWM, the DKI parameters showed better correlation to disability than DTI, suggesting that changes in kurtosis parameters may precede lesion formation. Standard DKI and DTI parameters produced results similar to those shown previously in other MS models. The estimated parameters of the WM model, however, yielded new microstructural information that could provide a key for improved understanding of EAE mechanisms.

Methods

Theory. Diffusion Kurtosis Imaging

DKI (Jensen et al., 2005) improves the approximation of the diffusion weighted signal in vivo (Filli et al., 2014; Raz et al., 2013; Rosenkrantz et al., 2015) and ex vivo (Veraart et al., 2011) by including the next term in the cumulant expansion (Kiselev, 2010; van Kampen, 2007) of the DWI signal S where D is the i,j element of the rank 2 symmetric diffusion tensor D and W is the i,j,k,l element of the symmetric rank 4 kurtosis tensor W, b is the diffusion weighting (b-value), and denotes the i-th component of measurement direction . In analogy to diffusion tensor based fractional anisotropy (FA), mean (MD), axial (Dǁ) and radial (D┴) diffusivity, the kurtosis tensor (W) provides additional biomarkers: kurtosis fractional anisotropy (KFA) (Hansen and Jespersen, 2016), mean of the kurtosis tensor (MKT) (Hansen et al., 2014, 2013) , axial (Kǁ) and radial 5 kurtosis (K┴) (Jensen and Helpern, 2010).

The choice of the DTI parameters assessed in this study was based on previous works. In particular, FA, MD, Dǁ and D┴ have been shown to be affected by MS pathology (Ceccarelli et al., 2007; de Kouchkovsky et al., 2016; Falangola et al., 2014; Guglielmetti et al., 2016; Inglese and Bester, 2010; Jelescu et al., 2016; Kelm et al., 2016; Mesaros et al., 2009). In DKI, a choice of MKT was motivated by a decrease in mean kurtosis that was shown in human GM and NAWM (Raz et al., 2013; Yoshida et al., 2013) and linked with cognitive impairment in MS (Bester et al., 2015). A decrease in axial kurtosis Kǁ and radial kurtosis K┴ was detected in an animal model of chronic MS (Falangola et al., 2014; Guglielmetti et al., 2016), K┴ was found to be related to the myelin content (Kelm et al., 2016).

In this work, we used a variant of the ‘standard’ WM model (WMM) that has been extensively explored recently (Fieremans et al., 2011; Jelescu et al., 2015a; Jespersen et al., 2007; Novikov et al., 2018b; Zhang et al., 2012). The model is designed to approximate diffusion inside and outside WM fascicles. It consists of two non-exchanging Gaussian compartments representing extra-axonal and intra-axonal space. Here, the diffusion in the extra-axonal space is described by an anisotropic cylindrically symmetric tensor which is defined by extra-axonal radial and axial diffusivities (De,┴ and De,ǁ). Axons, having radii much smaller than the diffusion distance, are assumed to appear as one-dimensional sticks and thus only the intra-axonal axial diffusivity D is non-vanishing. Taking f to be the volume fraction of the axonal compartment, and to be the fiber-orientation distribution function (fODF), the diffusion signal S measured in the direction can be written as:

In this study fiber bundles are not assumed to be parallel (as in (Fieremans et al., 2011)) but rotationally symmetric. Therefore, fODF follows the Watson distribution , where κ is the concentration parameter and is the symmetry axis. This model was chosen due to the fact that its axial symmetry assumption is valid for majority of spinal cord (SC) WM. In addition, it has a high range of validity, and is analytically related to DKI parameters (Jespersen et al., 2017; Novikov et al., 2018b).

For this study, the WMM parameters chosen to be assessed were those sensitive to neural damage (Falangola et al., 2014; Kelm et al., 2016); in particular, f, a biomarker for axonal loss (Fieremans et al., 2012) linked to myelin content and axon density (Kelm et al., 2016), D, which is associated with intra-axonal injury (Hui et al., 2012), De,┴ which through tortuosity is related to the g-ratio (Jelescu et al., 2016) and De,ǁ which is a marker of demyelination (Fieremans et al., 2012). In addition, a measure (Grussu et al., 2017) of fiber dispersion κ was studied.

For this study, DKI data was used as a starting point for WMM parameters estimation. This approach is commonly used in different models of tissue microstructure (Fieremans et al., 2011; Hansen et al., 2017; Hui et al., 2015; Jespersen et al., 2012; Novikov et al., 2018b; Novikov and Kiselev, 2010; Szczepankiewicz et al., 2016). Fitting DKI instead of straightforward fit of WMM parameters, enables usage of linear least squares algorithms that yield stable estimates (Chuhutin et al., 2017) decreasing the chances to end up in a local minimum. WMM fit at order b2 yields two solutions that fit data equally well (Novikov et al., 2018b). Using DKI fit and consequently estimating the WMM parameters allows to choose a particular solution branch (Jelescu et al., 2015b) explicitly.

Animal treatment

Female C57BL/6j bom (B6) mice aged 6–8 weeks obtained from Taconic Europe A/S, (Lille Skensved, Denmark) were maintained in the Biomedical Laboratory, University of Southern Denmark (Odense).

Mice were immunized at age of 8–10 weeks by injecting subcutaneously 100µl of an emulsion containing 100µg myelin oligodendrocyte glycoprotein (MOG)p35–55 (TAG Copenhagen A/S, Frederiksberg, Denmark) in incomplete Freund’s adjuvant (DIFCO, Alberstslund, Denmark) supplemented with 400 µg H37Ra Mycobacterium tuberculosis (DIFCO). Bordetella pertussis toxin (300 ng; Sigma-Aldrich, Brøndby, Denmark) in 200 µl of phosphate buffered saline (PBS) was injected intraperitoneally at day 0 and day 2. Animals were monitored daily from day 5 and scored as follows: 1, loss of tonus of up to half of the tail; 1.5, loss of tonus of more than half of the tail; 2, complete loss of tail tonus; 2.5, difficulty in walking (one leg); 3, Difficulty walking (both legs); 3.5, paresis in one hind leg; 4, paresis both hind legs; 4.5, paralysis of one hind leg; 5, paralysis of both hind legs. About 75% of the mice showed symptoms of EAE. All the scoring was performed by the same person (AW) with previous experience of EAE animal assessment (Wlodarczyk et al., 2014). Severe EAE usually developed 14–18 days after immunization. Based on the provided EAE-scale, the animals were divided into roughly equisized groups of samples: low-grade (EAE score 1.5–2, 5 samples), intermediate (2.5–4, 3 samples), high (4.5–5, 5 samples). If not stated otherwise, the control group (non-manipulated age matched animals) is henceforth referred to as zero-grade for convenience (5 samples).

Animal experiments were approved by Danish Animal Experiments Inspectorate (approval number 2014–15-0201–00369).

Sample preparation

Mice were euthanized by pentobarbytol overdose and transcardially perfused with PBS followed by 4% buffered paraformaldehyde (PFA) (pH 7.4). The spinal column was extracted and stored in 4% PFA for 7 days. On day 8, now fully fixed cords were manually dissected out of spinal column, and stored in 4% PFA with the meninges removed until MRI. 24 hours prior to the experiment, the samples were washed in PBS to remove PFA and to minimize associated T2*-related signal attenuation (Shepherd et al., 2009, 2005). The SC was cut into 3 parts and segments from T8 up to L6 were selected for imaging. We differentiate between three segments of mouse SC as follows: mid-thoracic (MTO):T8-T11, lower thoracic (LTO):T12-LU1, lumbar (LU):L2-L6.

MR imaging

Imaging was performed on a 16.4 T vertical bore Bruker Biospin (Ettlingen, Germany) Aeon Ascend magnet equipped with a Micro5 probe and a gradient unit capable of delivering up to 3000 mT/m in all directions. The samples were placed in a 5 mm NMR tube filled with Fluorinert© 3M and held parallel to the direction of the main magnetic field using a polypropylene straw. The temperature was monitored and maintained at 23.6°C using air flow.

In order to increase the prospective impact of the study we used a variant of standard Stejskal-Tanner pulse sequence that is commonly implemented in human scanners. Thus, diffusion kurtosis data were acquired using a 2D diffusion weighted fast spin echo sequence with echo train length (ETL=8), first echo time (effective echo time (TE) = 15 ms), echo spacing (ESP=4.23 ms), total repetition time (TR=2000ms)(Beaulieu et al., 1993; Kelm et al., 2016; West et al., 2018). Receiver bandwidth (BW) for signal acquisition = 83 kHz. For each SC, between 16 and 22 0.5mm-thick slices were scanned. For each slice, a matrix of size 120 × 120 voxels with field of view (FOV) (FOV=4.2 mm × 4.2 mm) (resolution of 0.035 mm × 0.035 mm) was acquired. Diffusion weighting was performed with short gradient pulses of duration (δ=1.5ms) and separation (diffusion time) Δ = 10ms. Diffusion weighting with b-values of 0.2,0.3,0.5,0.6,0.9,1,1.2,1.5,1.8,2.1,2.5 ms/µm2 was applied along 30 directions, with 1 average (NA) for b < 1.2 ms/µm2 and 2 averages for b > 1.2 ms/µm2. Sixty b = 0 ms/µm2 images were collected for normalization. The total scanning time per spinal cord was about 10 h. Examples of images acquired with b = 0.2ms/µm2 are provided in Fig. 1 (A,B). SNR (amplitude ratio) of the acquired raw data was estimated to be ~30–40 (in WM, b = 0), ~50 (in GM, b = 0), ~30 (in WM, b = bmax), ~20 (in GM, b = bmax).

Fig. 1.

Example of acquired raw data and a corresponding data fit. Subfigures (A) and (B) show an example of a raw signal image acquired for low diffusion weighting (b=0.2 ms/µm2) in mid thoracic and low thoracic segments of a control sample. Subfigure (C) shows a high grade sample, with a visible lesion acquired with the same low diffusion weighting (b=0.2 ms/µm2). Data and data fit that correspond to three different voxels in the slice denoted in (C) are shown in subplot (D). Lesion voxel location is marked in red, NAWM voxel in green and GM voxel in magenta in subplot (C). Multiple data points plotted under each b-value on the x-axis correspond to different directions. The inset shows the enlarged part of the graph corresponding to b=0.5 ms/µm2.

High resolution T2*-weighted images for lesion delineation were acquired using fast low angle shot (FLASH) pulse sequence with twice the in-plane resolution (0.018 mm × 0.018 mm) and the same slice thickness (0.5 mm), NA = 2 and TE = 5ms.

Image segmentation

Image segmentation of white and gray matter was performed manually based on the mouse spinal cord atlas (Watson, 2009).

Lesions were manually outlined on T2*-weighted slices as described in (Steinbrecher et al., 2005) and thereafter lesion contours were downsampled to the resolution of DWI maps. On each slice, potential abnormalities were inspected and compared to the atlas. Voxels with abnormal hyperintensity that could not be explained by the anatomical features of SC, were manually marked using an in-house developed software tool. Delineation followed a conservative definition of the lesion. As such, whenever there was a suspicion that the increase in WM intensity could be explained by anatomical features, the voxels were not delineated as lesion. The slices and spinal cords were presented in randomized order and the examiner (AC) was blinded to the grade. An example of this segmentation is shown in Fig. 2. NAWM was defined after the segmentation as a non-lesion WM. Lesion load was defined as fraction of volume taken by abnormal hyperintensity in T2*-weighted images.

Fig. 2.

Example outcome of lesion identification in four spinal cords in a lumbar segment. From left to right the grades are control, low grade, intermediate grade and high grade of EAE. For each of two subplots, the upper image represents a raw T2*-weighted image, while the lower image shows the same map with the manual lesion delineation superimposed in yellow.

Parameter estimation

The raw images were denoised using the Marchenko-Pastur PCA method (Veraart et al., 2015) and subsequently corrected for Gibbs ringing artefacts (Kellner et al., 2015) before further analyses. Standard single diffusion encoding technique that was used in this study is capable of accessing only fully symmetric diffusion and kurtosis tensors. Therefore, twenty-two independent components of fully symmetric diffusion and kurtosis tensors (Jensen et al., 2005) were fit to the data using interior-point weighted linear least squares algorithm (Veraart et al., 2013) implemented in Matlab. Based on (Chuhutin et al., 2017), WM voxels were fit up to a maximum b-value of bmax = 2.5 ms/µm2, and GM voxels were fit up to bmax = 1.2 ms/µm2. The fit quality was inspected for each sample. An example of data fit for a representative voxel in WM lesion, NAWM and GM is shown in Fig. 1 (C,D). Diffusion and kurtosis tensor parameters were calculated according to (Hansen et al., 2014, 2013; Jensen et al., 2005; Jensen and Helpern, 2010). The exact analytical derivations of WMM parameters from the elements of diffusion and kurtosis tensors used in this study are provided in (Jespersen et al., 2017) (assuming a Watson distribution of neurites). The general case is presented in (Novikov et al., 2018b). Different sets of WMM parameters can yield the same DKI signal, an effect known as parameter degeneracy, and a subject of current interest (Jelescu et al., 2015a, 2015b). However, in this work, only parameters corresponding to the so-called ‘plus’ branch (Hansen and Jespersen, 2017; Jespersen et al., 2017; Novikov et al., 2016), typically having D > De,ǁ were considered. This choice was guided by the recent publications (Kunz et al., 2018), that suggest ‘plus’ branch solution to be biologically plausible. Less than 10% of voxels in any slice displayed non-physical values, such as a negative diffusivity. These voxels were excluded from further statistical analysis of WMM parameters. In total, for all spinal cords, 245851 GM voxels and 246393 WM voxels were analyzed for DTI/DKI parameter estimation, while WMM parameters were estimated in 232274 voxels.

In WM, the estimated parameters were: axial diffusivity (Dǁ), radial diffusivity (D┴), fractional anisotropy (FA), axial kurtosis (Kǁ), radial kurtosis (K┴), which were included for a model independent assessment and the previously mentioned WMM parameters (extra-axonal radial De,┴ and axial De,ǁ diffusivities, intra-axonal diffusivity D, volume fraction of axonal compartment f, and concentration parameter of the Watson distribution, κ). In GM the low tissue anisotropy causes the estimated direction of primary eigenvector to be unstable/poorly defined, and thus, the values of axial and radial diffusivity and kurtosis are less reliable/meaningful. Due to that and in order to restrict the number of compared parameters to avoid unnecessary multiple comparisons, it was decided to limit the scope of estimated parameters in GM to MD and MKT.

Statistical Analysis

The voxels from all spinal cords were input to a linear mixed effects model (LME) (Gelman and Hill, 2007; Goldstein, 2011).

The choice of model was guided by current recommendations in (Barr et al., 2013; Bolker et al., 2009) and iterative maximization of Akaike information coefficient (Akaike, 1998).

The choice of LME to estimate and study the effects of EAE on kurtosis tensor parameters was guided by the observations that pathological changes both inside lesions, in NAWM and in GM contribute to clinical disability in both human MS and in animal models of neurodegeneration (de Kouchkovsky et al., 2016; Evangelou et al., 2000; Filippi et al., 2012; Filippi and Rocca, 2011; Inglese and Bester, 2010; Kipp et al., 2016; Lassmann and Bradl, 2016). The specific form of the LME designed to take into account random contributions of sample to sample variability. Our assessment of LME fitting quality (in Table 2) was in line with up-to-date recommendations for LME (Baayen et al., 2008; Edwards et al., 2008; Nakagawa et al., 2013).

Table 2:

The results of the fit of the linear mixed effects model. For each of the studied parameters (in rows), the following are presented in columns: percent of outlier values removed, quality of LME fit (Edwards et al., 2008) p-values for coefficients of grade, lesion, segment and segment:lesion, partial (Edwards et al., 2008) of the same 4 coefficients and the results of the FDR multiple comparison test. Since lesions were registered only in WM, the coefficients of lesion are absent in GM.

Tissue	name	Outliers %	Rβ2	P-values				partial Rβ,p2				Grade FDR
				Grade	Lesion	Segment	Segment :Grade	Grade	Lesion	Segment	Segment Grade
GM	MD	1.69	0.910	0.576		0.0598	0.382	0.117		0.588	0.295
GM	MKT	1.80	0.980	0.010		0.136	0.290	0.863		0.710	0.663	*
WM	K ┴	2.12	0.995	<0.001	<0.001	0.0111	0.004	0.953	0.991	0.896	0.839	*
WM	FA	3.50	0.977	0.042	<0.001	0.0023	0.004	0.852	0.982	0.972	0.852
WM	K ║	1.86	0.979	0.050	<0.001	<0.001	0.016	0.678	0.864	0.922	0.830
WM	D ┴	2.13	0.992	<0.001	<0.001	0.003	0.021	0.972	0.926	0.955	0.765	*
WM	D ║	1.72	0.991	0.590	<0.001	0.0034	0.378	0.132	0.903	0.838	0.408
WM	f	2.25	0.962	<0.001	<0.001	0.5822	0.011	0.967	0.971	0.259	0.853	*
WM	D a	1.52	0.971	0.525	0.123	0.0127	0.110	0.204	0.478	0.862	0.543
WM	D e,║	1.42	0.919	0.003	0.014	0.0044	0.349	0.886	0.738	0.923	0.412	*
WM	D e,┴	2.40	0.948	0.008	<0.001	0.0029	0.005	0.949	0.898	0.882	0.910	*
WM	κ	2.05	0.993	<0.001	<0.001	<0.001	<0.001	0.957	0.992	0.980	0.916	*

Each of 12 examined parameters p was thus fit to using Wilkinson notation (Wilkinson and Rogers, 1973), where g is grade, s is slice, lesion is l, a is sample (animal)[1]. The ‘fixed’ effects part of the model was designed to allow the parameters to depend on grade, while the size of the effect was permitted to be different in various SC segments (first term). The second term encodes the expected difference in parameter values inside (l = 0) and outside (l = 1) the T2* hyperintense lesions. Sample-to-sample variations were allowed by including ‘random’ effects for segment, grade and lesion, each grouped sample-wise.

To avoid a small number of data points having an undue influence on the regression, outliers 2.5 standard deviations above and below the model residual means were removed after the initial fit, and the model was refitted. This procedure is in agreement with literature (Baayen, 2008; Baayen et al., 2008; Tremblay and Tucker, 2011). By verifying that predictors are significant before and after the outlier removal we verified that the extreme values do not substantially influence the regression. The removed outliers constitute less than 4% of data. The comparison of in outliers and the rest of the data in terms of quality of fit (χ2) is given in Supplementary material.

For each of the ‘fixed’ effects, analysis of variance (ANOVA) p-values were calculated. These p-values represent the significance of individual fixed effects as well as the combined effect of segment and grade on parameter. The p-values describing the significance of the linear relationship between the measured parameter and the grade of disability of the EAE animal were finally reevaluated using the false discovery rate (FDR) procedure (Benjamini and Hochberg, 1995).

The quality of the fit of LME was estimated using (Edwards et al., 2008), so that where is a statistic corresponding to the null hypothesis H0: β1 = β2 = … = β = 0 for q – 1 fixed effects, β, ν is Satterthwaithe’s estimate of effective degrees of freedom. Partial were calculated to obtain the relative measure for each of the ‘fixed’ effects that corresponds to the null hypothesis H0 : β = 0 for j ϵ {1,…, q – 1}.

While the LME analysis was performed primarily to discern which DKI/WMM parameters correlate with disability and which tissues are affected, a following post-hoc analysis aimed to check the sensitivity of the chosen parameters. The parameter means were calculated for each sample, for each SC segment, and for the SC as a whole, in GM and WM separately. The parameters surviving FDR correction of LME p-values were checked and the grades with significantly different means were identified using two sample tests (one-way ANOVA, followed by FDR correction). For the post-hoc analysis, no outlier removal was performed.

A further post-hoc comparison of individual parameters inside lesions supported the initial LME model assumption of independence between the grade and lesion LME model parameters (data not provided). Thus, we based the post-hoc lesion analysis on the premise that the distribution of parameters inside lesions does not depend on segment or EAE grade.

Quantitative results (mean and standard deviation along the disability group, no outlier removal) for all the measured parameters are provided in Table 1. The most distinct changes are observed in GM MKT and NAWM K┴ (decrease) and an increase in De,ǁ in NAWM. Of all the measured parameters De,ǁ has the highest variance.

Table 1:

Fit results of all measured parameters averaged for disability group. For each parameter a mean estimate provided along with standard deviation of error inside the lesion (subplot A) and outside the lesion (in NAWM, subplot B). The mean and standard deviation of error was calculated only in the tissues in which a particular parameter was used for a successive LME analysis. Thus the statistics for MKT and MD was estimated only in GM, and the statistics for the rest of the DKI and WMM parameters was estimated only in WM. Note also that GM parameters (MD,MKT) were not calculated inside the lesions.

A:Lesions
Tissue	Parameter	Control		Low grade		Intermed grade		High grade
		Value	Stdev	Value	Stdev	Value	Stdev	Value	Stdev
GM	MD
GM	MKT
WM	K ┴	1.89	0.52	1.44	0.46	1.37	0.37	1.37	0.43
WM	FA	0.48	0.17	0.35	0.15	0.35	0.11	0.40	0.14
WM	K ║	1.04	0.17	0.96	0.14	0.89	0.12	0.94	0.12
WM	D ┴	0.47	0.20	0.52	0.17	0.52	0.12	0.45	0.10
WM	D ║	1.14	0.20	0.98	0.21	0.97	0.16	0.91	0.17
WM	f	0.54	0.08	0.51	0.09	0.49	0.07	0.48	0.08
WM	D a	1.93	0.34	1.98	0.38	1.98	0.27	1.86	0.26
WM	D e,║	0.45	0.18	0.50	0.25	0.55	0.18	0.49	0.17
WM	D e,┴	0.66	0.22	0.78	0.23	0.74	0.16	0.64	0.15
WM	κ	5.31	2.58	3.83	2.04	3.52	1.83	4.08	1.96
B:NAWM
Tissue	Parameter	Control		Low grade		Intermed grade		High grade
		Value	Stdev	Value	Stdev	Value	Stdev	Value	Stdev
GM	MD	0.57	0.06	0.58	0.06	0.58	0.06	0.58	0.07
GM	MKT	1.15	0.18	1.17	0.18	1.15	0.18	1.10	0.20
WM	K ┴	2.66	0.71	2.57	0.65	2.36	0.64	2.20	0.67
WM	FA	0.80	0.06	0.79	0.07	0.74	0.10	0.74	0.11
WM	K ║	0.91	0.14	0.87	0.13	0.89	0.14	0.87	0.15
WM	D ┴	0.21	0.05	0.23	0.06	0.26	0.08	0.26	0.08
WM	D ║	1.24	0.16	1.26	0.16	1.21	0.18	1.19	0.20
WM	f	0.64	0.07	0.64	0.06	0.61	0.06	0.61	0.07
WM	D a	1.94	0.22	1.95	0.21	1.90	0.21	1.91	0.24
WM	D e,║	0.36	0.19	0.38	0.19	0.39	0.16	0.41	0.19
WM	D e,┴	0.48	0.12	0.49	0.11	0.54	0.12	0.49	0.13
WM	κ	12.54	3.92	12.67	4.17	9.90	3.88	10.82	4.55

Results

Diffusion MRI: Parameter estimation

A one-way-ANOVA of sample-wise mean of relative lesion load showed that controls were significantly different from the diseased (EAE) animals in all segments (provided in Supplementary material). However, despite apparent increase in the mean values of lesion load, the variance of lesion load increases in SC slices of high grade animals. Subsequently, the lesion load is unable to statistically significantly discriminate between grades of disability.

Figure 3 shows maps of all the investigated parameters for a representative animal in each of the grades (control, low, intermediate, and high grade) in the medium thoracic (T9) segment. WMM parameters are restricted to the manually delineated WM to approximately fulfill the assumptions of the model. Qualitatively, the maps show an increase in asymmetry in animals with higher disability grade. The Watson concentration parameter κ displays the biggest variation in the maps. Parameters that are accessible with DTI (apart from MD) show better contrast between WM and GM. However, parameters derived from DKI and WMM show more variability inside WM.

Fig. 3.

A: Examples of parameter maps for each of the measured parameters in mid-thoracic segments of spinal cord. Each column (from left to right) corresponds to different grades of EAE disability: control animal, low grade, intermediate grade and high grade of EAE. Rows correspond to different measured parameters (from top to bottom): mean diffusivity, MKT, FA, axial diffusivity, radial diffusivity, radial kurtosis, parallel kurtosis

B Examples of parameter maps for each of the measured parameters in mid-thoracic segments of spinal cord. Each column (from left to right) corresponds to different grades of EAE disability: control animal, low grade, intermediate grade and high grade of EAE. Rows correspond to different measured parameters: axonal water fraction, axonal diffusivity, axial extra-axonal diffusivity, radial extra-axonal diffusivity and concentration parameter of Watson distribution, the upper row depicts the delineation of spinal cord on the background of FA map.

Diffusion MRI: validating metrics with LME

Table 2 shows the estimators of the LME fit quality and values that quantify the capability of LME model parameters to explain each of the parameters.

All parameters demonstrated relatively good quality of fit . MD and De, attained the lowest values of ~0.91.

P-values shown in Table 2 quantify the extent to which each of the 12 studied parameters can be explained by the parameters of the linear mixed effects model. Grade had a significant effect on 7 out of 12 parameters after FDR: MKT in GM, and K┴, D┴, De,, De,┴, κ and f in WM. All parameters except MD and MKT in GM and f in WM were found to depend significantly on the segment. Likewise, the interaction between grade and segment was statistically significant in all but five parameters, i.e. the two GM parameters MD and MKT, three WM parameters De,, D and D. All parameters but D were significantly different between lesion and normal appearing brain tissue.

The results of the calculation of partial (Edwards et al., 2008) for each fixed effect variable are also provided in Table 2. revealed high association between the kurtosis/WMM parameters and the disability grade of EAE (~0.9) for all the parameters that were found significant in FDR procedure. The comparison of values showed that the disability grade accounted for most of the variation in 4 out of 12 parameters D┴, De,ǁ, De,┴ in WM and MKT in GM.

Additional characteristics of the LME fit are provided as Supplementary material. These include a different measure of fit quality (Johnson, 2014; Nakagawa et al., 2013) and estimates and confidence intervals of fixed effects of LME-model. These estimates show that among the parameters which are significantly correlated with the grade, MKT, K┴, f, De,┴ and κ decrease with the increase in disability grade, while D┴ and De,ǁ increase with increasing grade. The LME results have the same direction and magnitude of decrease or increase as overall means in Table 1.

Diffusion MRI: Post-hoc statistical analysis

From the LME analysis, we found that the variation of several GM and WM parameters can be explained by EAE-grade and by lesion status, i.e. whether or not the voxel is located inside a lesion. A follow-up post-hoc analysis intended to investigate group-wise behavior of the segment-wise means in parameters with a significant grade. In particular, Table 3 shows the results of the post-hoc analysis of sample means outside the lesions, in GM and NAWM.

Table 3:

Post-hoc analysis of parameters in GM and NAWM. Average value for all the NAWM or GM in the particular segment in each sample. (A) For each of the disability groups comparisons low-grade vs intermediate grade, low-grade vs high grade and intermediate grade vs high grade parameters that were found significant (p<0.05) after ANOVA of per-sample mean in each one of the segments is provided in the corresponding cell. An FDR correction of multiple comparisons has been taken into account. GM parameters are underlined. (B) Each of disability groups (low, intermediate, high) compared with the control group. Parameters that were found significant (p<0.05) after ANOVA of per-sample mean in each of the segments are listed in corresponding cells. An FDR correction has been performed. GM parameters are underlined.

A
	Low grade	Intermediate grade	High grade
Low grade		MTO: MKT	MTO: MKTK┴ f De,║LTO: K┴ fLU: K┴ D┴ fAll: K┴ f
Intermediate grade			MTO: K┴
B
	Low grade	Intermediate grade	High grade
Control	MTO:	MTO:	MTO: MKTK┴ f De,║
	LTO:	LTO:	LTO: K┴ f De,║
	LU:	LU: FA K┴ D┴ f	LU: FA K┴ D┴ f De,║ κ
	All:	All: De,║	All: K┴ D┴ f De,║

In GM, MKT showed significant difference between the grades only in the slices located in mid-thoracic segment of the spinal cord, specifically between low and intermediate and between low and high grades.

In NAWM, 5 out of the 6 biomarkers surviving FDR correction demonstrated significant differences between the control and diseased animals, mainly in the lumbar SC. Two DKI parameters (K┴, D┴), and two WMM parameters (f, De,ǁ) depended significantly on EAE grade. K┴ and f were found to survive pooling all segments together, demonstrating an overall significant difference between high and low grade.

As an illustration, a representative part of the data and the associated post-hoc analysis is given in Fig. 4. In this figure, the distributions over voxels of MKT in GM (Fig. 4 (A)) and K┴ in NAWM (Fig. 4 (B)) are visualized using box plots for each of the grades and segments. The means of each of the spinal cords, which were used in the post-hoc analysis (Table 3), are superimposed on the boxplots as blue circles. The significantly different grades are marked by asterisks. Note that the seemingly large number of outliers apparent in the box-plots of Fig. 4 constitute a small fraction of the more than 10000 voxels sampled for each spinal cord.

Fig. 4.

Examples of parameter distributions in the post-hoc parameter analysis for MKT in GM (A) and K┴ in NAWM (B) illustrated with box-plots. Each of 4 subplots corresponds to one of the spinal cord segments (mid-thoracic, low-thoracic and lumbar and a graph for all voxels pooled across segments). In each one of the subfigures a box-plot represents the distribution of values inside and outside the hyperintensity lesions. The outliers (voxels) are plotted individually as red crosses. Blue dots correspond to the parameter means within each spinal cord. Asterisk denotes significant group-wise difference between the spinal cord means. In all three plots the central mark indicates the data median, the bottom and top edges indicate the 25th and 75th percentiles. The whiskers extend to the most extreme data points excluding outliers.

The same type of post-hoc analysis that was used to study the voxels outside of the lesions, was used to investigate voxels inside the lesions. The analysis revealed that the vast majority of segment-wise means inside the hyperintensity lesions did not show any significant differences between EAE-grades, with only Kǁ in lower thoracic segment showing difference between the grades at the level of significance p < 0.05 (results of this analysis are provided as Supplementary material).

The results of an additional post-hoc analysis adapted to test for difference in parameter means between lesions and NAWM, revealed that the difference is significant for all 10 WM parameters are provided in Supplementary material.

Combining T2*-weighted images and diffusion MRI

Based on the presented results, that report lesion load to be insignificantly related to disability grade, we next consider a “hybrid” way of addressing the lesion load and DKI parameters. It can be done e.g. using a compound variable that reflects both lesion load and NAWM health to provide a way of distinguishing different grades of EAE. In particular, an animal-wise LME fit of the model where g is the grade and l is lesion load, and are animal-wise mean values of De,ǁ and K┴ in NAWM, yielded coefficient values p0 = 2.9, p1 = 12.7, p2 = −2.3, p3 = 10.4 ms/µm2. A “hybrid” metric that uses these parameters is able to distinguish not just between control and high, control and intermediate groups but also between low and high, low and intermediate EAE-grades. The results of using such metric are shown in Fig. 5. However, a follow up study that will test this hybrid metric with an independent data is needed.

Fig. 5.

Application of a “hybrid” biomarker (Eq. 6) on animal-wise data. From left to right the values are group-wise means of described control, low-grade, intermediate and high grade EAE. Error bars depict standard deviation of animal-wise estimates of the biomarker. Asterisk denotes statistical significance measured with 1-way ANOVA.

Discussion

Mapping quantitative biomarkers in MS – as well as in other neurological disorders – is essential for early diagnosis, follow up of treatments, and testing novel avenues for treating disease. Lesion load is a biomarker that is easily and commonly imaged both in humans and in preclinical studies. However, even though it has been historically associated with motor deficits in MS and EAE (Bjartmar et al., 2001; Sathornsumetee et al., 2000), its correlation with disability is poor (Barkhof, 2002; Bergers et al., 2002; Robinson et al., 2010; Wuerfel et al., 2007), a disparity known as the clinico-radiological paradox (Cohen et al., 2016; Nathoo et al., 2014). Confirming the clinico-radiological paradox in MS and previously published results in EAE (Wuerfel et al., 2007), we did not find correlations between the EAE-grades and lesion load. In addition, using DKI provided no obvious advantages for structural segmentation. We therefore turned to LME fitting of the metrics that map microstructural aspects of the tissue. Instead of trying to improve on lesion load detection, the study aimed to test additional quantitative metrics obtained from DKI against disease severity in an animal model of MS.

In GM, MKT depended significantly on disability grade (Tables 2 and 3; Fig. 4). This is in line with human studies (Agosta et al., 2007; Bester et al., 2015; Raz et al., 2013; Zackowski et al., 2009) reporting similar changes in GM. Such changes are likely indications of GM pathology, possibly associated with neuronal degeneration and myelin loss (Guglielmetti et al., 2016). Interestingly, while the biggest burden of lesions and most of the changes in NAWM (Table 3; Fig. 4) were associated with the lumbar section of the spinal cord, most of the changes detected in GM were observed in mid-thoracic sections. Given that the EAE-induced disability progresses from hind- to forelimbs, one might hypothesize that the damaged GM tissue in mid thoracic SC is connected to the damaged fascicles in lumbar WM. Thus, in EAE, correlated pathological mechanisms may be responsible for damage in NAWM and GM. This may be similar to human MS, where the spatial and temporal relationships between the damage in GM and NAWM are still not fully resolved and might depend on disease phenotype (Bodini et al., 2009; Pirko et al., 2007; Steenwijk et al., 2015; Tewarie et al., 2018).

In NAWM, K┴ and D┴ showed the strongest and most robust results among tissue biomarkers derived from kurtosis and diffusion tensors. Radial kurtosis showed the strongest inverse relationship with EAE grade (i.e. it decreased with increasing disease severity). Such changes have also been observed previously in preclinical models of MS (Falangola et al., 2014; Jelescu et al., 2016; Kelm et al., 2016), while the opposite effect was observed in (Guglielmetti et al., 2016). This agreement might suggest closer similarities of the EAE to cuprizone or genetically induced chronic demyelination than to the acute inflammatory demyelination used in (Guglielmetti et al., 2016). An increase in D┴ was also found to be significantly correlated with EAE grade. Again, this behavior agrees with previous chronic demyelination studies (Falangola et al., 2014; Jelescu et al., 2016; Kelm et al., 2016) and with results of a previous DTI EAE study (Budde et al., 2009). Early human studies demonstrated similar behavior of D┴ and associated it with demyelination (Klawiter et al., 2011) and possible axonal loss (Naismith et al., 2010).

Among WM model parameters, De,ǁ was the one affected the most by EAE grade, while De,┴ was affected in a much weaker manner and with no significant effects in post-hoc analysis (Tables 2 and 3). Counter-intuitively, an increase in De,ǁ with grade was found. This fact could potentially be explained by axonal damage, changes in the structure of glial cells, and myelin loss, causing the extra-axonal space to have higher diffusivity. This result is in contrast with cuprizone models (Falangola et al., 2014; Guglielmetti et al., 2016; Jelescu et al., 2015b). The disparity can stem from differences between the mechanisms underlying tissue degeneration but also from other microstructural differences between the neural tissue in cerebrum and in spinal cord. Alternatively, it may also be a result of choosing a different solution ‘branch’ when finding parameters of WMM model (Hansen et al., 2017; Jelescu et al., 2015b; Jespersen et al., 2017).

The axonal water fraction (previously suggested to be a biomarker of axonal loss (Fieremans et al., 2011)) was also significantly affected by the differences in EAE grade. We found f to decrease with an increase in EAE-grade, therefore an axonal loss in NAWM could be one driver of the disability.

The ratio λ = De,ǁ/De,┴ (tortuosity) has been proposed as a biomarker of demyelination (Fieremans et al., 2012). In our data (see Results and Supplementary material), De,ǁ increased strongly and De,┴ decreased, thus overall tortuosity increased with increase in grade (the effect of grade on tortuosity variation was found to be significant in post-hoc LME fit, data available upon request). This finding is in contrast to (de Kouchkovsky et al., 2016; Falangola et al., 2014) and may provide evidence of pathological processes in EAE.

Our data shows that the Watson concentration parameter κ significantly decreased in a way that could be explained by EAE grade. This might be a result of axonal damage that could cause the breaking of fascicles and fanning out of individual axons. According to post-hoc analysis, this behavior was present in the lumbar segment of the spinal cord. Our result is in line with (Schneider et al., 2017), where increased fiber dispersion in NAWM of MS SC was reported. However, in our study increased fiber dispersion is present inside the lesions in mouse SC as well. This is in contrast with (Grussu et al., 2017), where a decrease in neurite orientation dispersion was measured in lesions of MS using neurite orientation dispersion and density imaging (NODDI) and histology. One possible reason for the disparity is different species (animal and human), where different pathological mechanisms could be at play. Another reason could be related to focal lesions being not as well defined in rodent models, and in EAE in particular, as in humans. Therefore, even though lesion detection was performed in a consistent and ‘blind’ way, a systematic error may have been introduced, e.g. if too big portions of NAWM are segmented as lesions. Third, differences in the employed diffusion models could be responsible for the disparity. This discrepancy may be addressed by validation of fiber dispersion in EAE and MS lesions using microscopy.

Dǁ has previously been shown to decrease significantly with EAE score and with axonal injury (Budde et al., 2009). Both Kǁ and Dǁ were significantly affected in some cuprizone studies (Falangola et al., 2014; Guglielmetti et al., 2016) but not in (Jelescu et al., 2016; Kelm et al., 2016). Our study found no evidence of any correlation between EAE grade and Dǁ or Kǁ. Consequently, our work provides an indication that in EAE, tissue changes due to demyelination and axonal loss are insufficient to change diffusivity or kurtosis parallel to fiber bundles.

In our study, FA showed no significant dependence on grade. This observation is in line with (Guglielmetti et al., 2016) where FA was not able to differentiate between the treatment groups and control.

Both the results of LME model fitting and post-hoc analysis demonstrated that all parameters but one (D) were able to distinguish lesions from NAWM. Together with that, in a more intriguing finding of this work the vast majority of the biomarkers measured within lesions did not differ across grades. A probable explanation of that is that there is no detectable differences in lesion microstructure between the grades. This may be evidence that CNS tissue that makes up WM lesions does not contribute to disability in EAE.

Recent studies (By et al., 2018, 2017; Grussu et al., 2015; Schneider et al., 2017) have applied NODDI and spherical means (SMT) techniques to spinal cord tissue WM in healthy controls and in MS patients and demonstrated promising diagnostic results. However, our findings suggest that the basic assumptions of these studies may not be fulfilled in all cases. In particular, since D was found not to be driven by EAE grade while De,ǁ increased with grade, those two parameters cannot be constant as assumed in NODDI (D = De,ǁ = 1.7µm2/ms). The constant tortuosity assumption (imposed in both NODDI and SMT) was also not valid in our dataset. Thus, assuming that our results are valid in human MS, the values estimated with those two techniques would be unable to reveal the true microstructural changes associated with disease progression and disability. Thus, it may prove necessary to allow all the parameters of the so-called ‘standard’ model (Novikov et al., 2018b; Jelescu et al., 2015a; Jespersen et al., 2007) to be determined by the fit.

Limitations

In this work, fixed tissue was used, which allowed longer scanning and better data quality (e.g due to motion and susceptibility effects of vertebrae) in comparison with in vivo protocols. This choice was justified by the assumption that despite known impact of fixation on tissue properties (Shepherd et al., 2009, 2005; Sun et al., 2005) the pathological effect on damaged tissue will be strong enough to be detectable in the present exploration study. Nonetheless, the results of this study are influenced by differences between ex-vivo and in-vivo tissues, which have not been fully accounted for yet (Horowitz et al., 2015).

Manual lesion segmentation, even though it was performed ‘blindly’, can potentially result in a systematic bias in contrast between estimated NAWM and lesion values, as well as increased variability.

Concerns have been expressed (Lampinen et al., 2018) about applicability of modeling constraints required for multi-compartmental modeling of diffusion in neural tissue. The compartment models for diffusion used in this work have not been fully validated across different tissue types, in vivo and ex vivo datasets, etc. However, in this study the modeling efforts were restricted to ex-vivo mouse spinal cord, where, according to histology (Ong et al., 2008) the axonal size is around 1µm and in the chosen regime of gradient strengths and waveforms the attenuation due to diffusion along the diameter is negligible (Dyrby et al., 2013; Nilsson et al., 2017). Chosen diffusion times are short enough (10ms), so that the exchange has only a minimal effect (Nilsson et al., 2013, 2009). Therefore, there is good reason to believe that in this case, the approximation of axons as sticks is approximately valid. Further investigation of the validity of the attained results e.g. the role of myelin water and related compartment-dependent T2-relaxation of diffusion weighted signal is necessary before translating the results of this study to the clinic.

Clinical implications

This work shows that the disability in EAE and therefore probably also disability in MS is correlated with and maybe is even driven by the neural matter outside the lesions.

In line with other works, this study also suggests the prominent role of SC GM (Agosta et al., 2007; Bester et al., 2015; Guglielmetti et al., 2016; Raz et al., 2013; Zackowski et al., 2009) and NAWM (Falangola et al., 2014; Jelescu et al., 2016; Kelm et al., 2016) in the development of disability. All these results support the fact that the search for disability biomarkers in human SC should concentrate on the neural matter outside the WM lesions. Future longitudinal studies could elucidate whether GM damage is a precursor of damage in NAWM.

This study was performed with 11 b-values and 30 directions, yielding approximately 350 images. However, such a big number of b-values was primarily needed to estimate parameters in different tissues. If a similar approach would be used but focussed e.g. only on the WM, the same type of analysis can be performed with 5 b-values and 30 directions, which can be achieved within a clinically feasible scan time of around 15 minutes. At the same time, due to relatively low b-values, this protocol is more accessible for clinical systems than two compartment models such as CHARMED (Assaf et al., 2004; Assaf and Basser, 2005; Barazany et al., 2009; De Santis et al., 2016) developed to attempt axonal diameter mapping.

This study points on the perspective of using WMM, where extracellular parallel diffusivity and axonal water fraction are recommended for assessment of human MS disability using SC MRI. Parameters of the full WMM in human MS spinal cord can be estimated and compared to our results along with other animal models (e.g. progressive model as in (Al-Izki et al., 2012)). Comparing between the biomarkers in different models can help to reassess models of spinal cord pathology in MS. Since the damage in spinal cord is better correlated with accrual of long-term disability (Inglese and Bester, 2010; Lin et al., 2006) than damage in cerebrum, the results of such an assessment can improve the understanding of mechanisms of MS progression.

There are some barriers in translation of the results of this work to a clinical setting. In addition to human MS pathology being distinct from the EAE animal model in terms of illness onset and its progression, human scanners feature multiple technical differences compared to the system used here. Such differences pose some challenges in the adaptation of the described methods.

In particular, this study has been performed with an ultra-high magnetic field strength. On clinical systems with much lower magnetic field, the relaxation times are different. Since myelination and compartmental relaxation significantly influence the parameters of microstructural model (Lampinen et al., 2019), the magnitude of the effect of disability on WMM parameters may be altered. In addition, robust determination of model parameters, is not straightforward on clinical scanners in vivo, mainly due to gradient strength limitations, and this is an active area of research. Current promising paths involve acquisition of “orthogonal information” that constrains parameter estimation. Such information can be provided by unconventional diffusion sequences e.g. Double Diffusion Encoding (Coelho et al., 2019; Reisert et al., 2019). Until the effect size of disability on our white matter model has been investigated on clinical systems, we recommend to use the DKI-parameters: K┴, D┴ in NAWM and MKT in GM.

An additional obstacle in translation of the presented results to clinical setting is due to motion artefacts and susceptibility effects of vertebrae. The state-of-the-art method (Wheeler-Kingshott et al., 2014) to mitigate these is using cardiac triggering and careful shimming by e.g dedicated spinal cord shimming coils (Topfer et al., 2016). Recent studies (By et al., 2018, 2017; Hori et al., 2014; Li and Wang, 2017; Raz et al., 2013), which were designed based on these recommendations, succeeded to estimate DKI biomarkers in human cervical spinal cord. In addition to that, dedicated acquisition methods (Hansen et al., 2016) can speed up the acquisition of DKI metrics (such as MKT in GM and K┴ in NAWM) that were found sensitive in this study.

Overall the DKI data as acquired in human scanners has lower SNR and is more artefact-prone than in preclinical ex vivo setup, and therefore there are still considerable issues to be solved before applying these methods on the MS patients. However, the high significance of changes presented in this study suggests that with adoption of new methods, these results will be relevant for clinical use.

In NAWM and GM the relationship between the disability and DKI and DTI metrics was found to be similar to other hypomyelinating MS models and to ex-vivo MS tissue.

In NAWM, changes in WM-modeling parameters (strong increase in De,ǁ, weak effect in D, De,┴) were clearly different to what has been observed in other animal models of MS.

The statistical analysis based on linear mixed effect models disentangled NAWM and lesion effects. Neither accumulated lesion load, nor DWI biomarkers in the tissue restricted by T2*-weighted lesion show any significant effect of lesions on EAE-grades.

A strong increase in De,ǁ of NAWM is an effect that has not been previously observed in other models of MS.

Summary

Overall, this study detects significant alterations in NAWM and GM (but not in WM lesions) in SC of EAE animals that correlate with dysfunction. These alterations are best detected with DKI and WMM biomarkers.