A three year follow-up study of gadolinium enhanced and non-enhanced regions in multiple sclerosis lesions using a multi-compartment T2 relaxometry model

Demyelination, axonal damage and inflammation are critical indicators of the onset and progress of neurodegenerative diseases such as multiple sclerosis (MS) in patients. Due to physical limitations of imaging such as acquisition time and imaging resolution, a voxel in a MR image is heterogeneous in terms of tissue microstructure such as myelin, axons, intra and extra cellular fluids and free water. We present a multi-compartment tissue model which estimates the water fraction (WF) of tissues with short, medium and high T2 relaxation times in a T2 relaxometry MRI voxel. The proposed method is validated on test-retest data of healthy controls. This model was then used to study longitudinal trends of the tissue microstructures for two sub-regions of the lesions: gadolinium enhanced (E+) and non-enhanced (L–) regions of MS lesions in 10 MS patients over a period of three years. The water fraction values in E+ and L– regions were found to be significantly different (p < 0.05) over the period of first three months. The results of this study also showed that the estimates of the proposed T2 relaxometry model on brain tissue microstructures have potential to distinguish between regions undergoing active blood brain barrier breakdown from the other regions of the lesion.


Introduction
Magnetic resonance imaging (MRI) is one of the most widely used in-vivo imaging 2 method for obtaining information on brain health. However, MRI voxels have limited 3 resolution due to physical constraints. Use of advanced MRI techniques such as T 2 4 relaxometry and diffusion weighted MRI can address this issue by providing 5 quantitative estimates on brain tissue microstructures such as myelinated and 6 unmyelinated axons, brain fiber orientation etc [1]. In patients with neurodegenerative 7 diseases (such as multiple sclerosis (MS)), obtaining information on health of axons and 8 myelin using in-vivo imaging techniques can improve our understanding of the current 9 status and the progress of the disease in the patients. Obtaining this information can 10 provide critical knowledge on the brain health. 11 The myelin, axons and fluids in the brain can be distinguished based on the nature 12 of their water content. Myelin is a tightly wrapped structure around the axons. The 13 water in myelin is very tightly bound to its surface. Hence myelin has a very short T 2 14 relaxation time compared to the other brain tissue structures. The axons, gray matter 15 cells and other cells in the brain have a higher T 2 relaxation time than that of myelin 16 but less than that of the fluids. Free fluids have the largest T 2 relaxation time. Hence, 17 based on the T 2 relaxation time brain tissues can be broadly classified into three 18 categories, short−, medium− and high−T 2 components [2,3]. Since a voxel in a MR 19 image of the brain is heterogeneous, there is a certain amount of the three T 2 20 components present in each voxel. In addition to these, in presence of an infection there 21 might also be some fluid accumulation. Hence a quantitative metric that conveys 22 information on the condition of these tissues can provide useful insights into the current 23 brain tissue health. 24 Quantitative T 2 relaxometry MRI sequences offer the capability to distinguish 25 tissues based on their T 2 relaxation times. T 2 relaxometry MRI has been used 26 effectively to calculate the myelin water fraction using a variety of approaches. Whittall 27 et al. and MacKay et al. [2,4] obtained myelin water fraction using a multi-component 28 T 2 model [5]. The multicomponent T 2 model describes the observed T 2 decay curve as a 29 weighted sum of an arbitrary number of decay curves. The weight of each T 2 decay 30 curve is obtained using non-negative least squares (NNLS). Whittall et al. obtained 31 myelin water fraction of the brain tissues by assuming the T 2 values of water in myelin 32 to be in the range of 10 − 55 milliseconds [2]. The multi-component T 2 models select an 33 arbitrary number of T 2 decay curves to model the observed decay signal [5]. For 34 example, Whittall et al. [2] considered 80 logarithmically spaced values between 15 35 milliseconds and 2 seconds. Since the number of parameters to be estimated 36 considerably outweighs the number of observations (i.e. number of echoes), some 37 regularization is mandatory. However, the choice of the regularization term and its 38 extent affects the fitting measure [6]. 39 An alternate approach is a multi-compartment T 2 relaxometry model where the 40 different T 2 water pools in brain tissues are modeled as a weighted mixture of 41 pre-defined continuous functions. Stanisz and Henkelman [7] fitted two T 2 components 42 (short and long T 2 in bovine optic nerve), each having a gaussian distribution on a 43 logarithmic scale. Akhondi-Asl et al. proposed a framework to calculate myelin water 44 fraction by modeling the R 2 space as a weighted mixture of inverse Gaussian 45 distribution of fast−, medium− and slow− decaying components with respect to T 2 46 relaxation times [8]. Layton et al. proposed a maximum likelihood estimation approach 47 to estimate the myelin water fraction values [9]. They obtained the Cramér-Rao lower 48 bounds of the variables to establish the difficult nature of simultaneously estimating the 49 model parameters and weights for each compartment for such models. In contrast to the 50 multiple exponential T 2 fitting approach, these models do not suffer from the ill-posed 51 weights estimation problem. In these models, the regularization is performed a priori 52 rather than a posteriori.

53
In this work, we propose a method for computing water fractions corresponding to 54 fast−, medium− and slow− decaying components with respect to T 2 relaxation times 55 from T 2 relaxometry MRI data acquired using 2D multislice Carr-Purcell-Meiboom-Gill 56 (CPMG) sequence. In this estimation framework, the T 2 space is modeled as a weighted 57 mixture of three continuous probability density functions (PDF) representing the three 58 T 2 compartments. For such models, the robustness and accuracy of the implementations 59 to simultaneously estimate the weights and parameters of the distributions was found to 60 be non-trivial and not reliable [9]. We therefore chose to take advantage of the earlier 61 studies performed in this field to fix the parameters of the distributions representing the 62 T 2 relaxation components [3,10,11]. Hence only the estimation of weights remain, 63 making the estimation model robust. In our work, we need to correct for the effect of 64 the stimulated echoes due to imperfect refocusing (due to the B 1 inhomogeneities) that 65 leads to errors in T 2 estimation [12]. The EPG algorithm is used to account for these stimulated echoes [13]. The field inhomogeneity (B 1 ) is estimated numerically. Since

72
In the last part of this work, we show an application of the method on MS patients. 73 We report observations on the evolution of the quantitative markers obtained from the 74 model in enhancing (MS lesion regions undergoing active blood brain barrier breakdown) 75 and non-enhancing MS lesion regions. Vavasour et al. studied the MWF and total water 76 content in three new MS lesions in two patients over a year [14]. Levesque et al. studied 77 evolution of MWF and geometric mean of T 2 values of gadolinium enhancing lesions in 78 five MS patients over a period of one year [15]. Although MWF was found to be quite MWF is a combined measure of edema and demyelination [16]. The MWF is a relative 86 measure and a change in its values should be studied in conjunction with the remaining 87 water fraction measures. In our analysis, we observed and compared the change in water 88 fractions of tissues with short, medium and high T 2 values in enhancing and with parameters p. As mentioned in the introduction, for robustness reasons, the PDF 100 parameters are pre-selected and fixed keeping in mind histology findings reported in the 101 literature [3,11]. In Eq. (1), M 0 is the magnetization constant. EP G(·) represents the 102 stimulated echo computed at the time point (t i = i × T E) using the EPG the voxel at time t i is expressed as: The parameters to be estimated in Eq. (2) are: {α j } 3 j=1 , B 1 . The optimization is thus formulated as a least square problem: where m is the number of echoes; Y ∈ R m is the observed signal; The parameters to be estimated in the least squares optimization problem stated in 112 Eq. (3), α and B 1 , are linear and non-linear in nature respectively. However they are 113 linearly separable. Hence optimization for α and B 1 is performed alternatively until 114 convergence is obtained in a desired error limit. In the first step, α is computed by 115 non-negative least squares (NNLS) optimization [17] with a fixed B 1 value. In the next 116 step, the weights computed in the first step are used to compute B 1 by a gradient free 117 optimizer (BOBYQA). We choose to perform a numerical optimization to obtain B 1 as 118 it does not have any closed form solution [13]. The integral in Eq. (4) also does not

129
Repeatability is an important aspect of quantitative MRI techniques. In the first 130 experiment we carried out test-retest experiments to assess the repeatability of the 131 proposed method. An important motivation of this study is to gain more insights into 132 the neurodegenerative disease phenomenon using the tissue microstructure information. 133 Hence, in the second experiment we performed a 3-year follow-up study on 10 multiple 134 sclerosis patients.  Healthy controls The age of the healthy controls was in the range of 26-32 years. 15 141 regions of interest (ROI) were marked in the brain for each healthy control over which 142 the test and retest values of the compartments' water fractions were compared. An 143 illustration of these ROI on a subject is shown in Fig. 1.
Assessing repeatability from test-retest data Test retest scans were performed for 4 healthy controls to study the repeatability of the proposed method. This figure shows the 15 regions which were marked on the healthy controls over which the repeatability was studied.   Multiple sclerosis lesions are focal lesions and grow in a concentric manner [18]. In the 155 early stages, brain tissues in the MS lesions undergo active blood brain barrier 156 breakdown [18,19]. Surrounding the core of the lesion is the edema as a result of tissue 157 inflammation due to ongoing tissue damage. As compared to the normal appearing 158 brain matter, the entire MS lesion regions appear hyper-intense on T 2 weighted MRI.

159
However, only the regions of the lesion undergoing active blood brain barrier breakdown 160 appear hyperintense on T 1 weighted MRI acquired post Gd injection. Hence, lesions in 161 active state have two regions, a region which actively undergoes blood brain barrier 162 breakdown and the regions which do not.

163
The objective of this experiment is to study the evolution of compartments' water  The lesion ROIs are marked at baseline and the same region is observed over a period of 193 three years. In the 10 MS patients, we observed 229 L− and 25 E+ lesion regions.

194
Since the lesions were marked on T 2 weighted images, all processed images were 195 registered to the T 2 images using a block matching algorithm [20,21]. The Bland-Altman (BA) plots for short, medium and high-T 2 water fraction estimates 199 over the ROIs are shown in Fig. 3. BA plots are scatter plots between the average   Table 1. observed. The gray area around the mean error level is its 95% confidence interval (CI  Fig. 3 for short, medium and high-T 2 water fraction estimates.

216
In this experiment we observed and compared the evolution of water fraction maps of  An example illustrating the comparison between the water fraction maps for a 225 healthy control and MS patient is shown in Fig. 4. Lesion-1 in the MS patient has a 226 large active core, whereas a very small region of lesion-2 is active. Both lesions show 227 indications of extensive demyelination. The medium-T 2 and high-T 2 water fraction 228 maps show varying trends among the lesions, and also between the regions of the lesion 229 undergoing active blood brain barrier breakdown and otherwise.

230
Short T 2 water fraction Results for short-T 2 water fraction (w s ) maps is shown in 231 Fig. 5. The w s values of the E+ and L− lesions at all time points is shown in Fig. 5(a). 232 The L− lesion regions are significantly associated with higher w s values as compared to 233 E+ at M00 (p = 0.014). However, the w s distribution of E+ and L− regions at the end 234 of 3 years are similar. The results of the change in w s values (∆w s ) between consecutive 235 scans in shown in Fig. 5(b). Largely negative ∆w s,0 values for E+ lesion regions  Fig. 5(a). The change in w s between consecutive scans is shown in Fig. 5(b). Significant group differences between groups are shown using * (p < 0.05) and * * (p < 0.01).

240
Medium-T 2 water fraction The evolution of the medium-T 2 water fraction (w m ) in 241 E+ and L− lesion regions is shown in Fig. 6(a). Although the w m values for both 242 groups reduce slightly at the end of 3 years, there is no evidence of difference between 243 E+ and L− with respect to the change in w m values between successive scans (refer 244 Fig. 6(b)).

245
High-T 2 water fraction The high-T 2 water fraction (w h ) values for E+ and L− 246 lesion regions are shown in Fig. 7(a). The E+ lesion regions are significantly associated 247 with a higher value of w h as compared to the L− population at M00 (p = 0.002).

248
Largely negative ∆w h,0 values (refer Fig. 7(b)) for E+ lesion regions suggest a decrease 249 in their w h values between M00 and M03. Subsequently, E+ lesion regions undergo  Fig. 6(a). The change in w m between consecutive scans is shown in Fig. 6(b). Significant group differences between groups are shown using * (p < 0.05) and * * (p < 0.01).  Fig. 7(a). The change in w h between consecutive scans is shown in Fig. 7(b). Significant group differences between groups are shown using * (p < 0.05) and * * (p < 0.01). 0.011 34.55% The common language (CL) effect size are reported for the measurements which were significantly different for E+ and L−. In this work, the CL effect size denotes the percentage of times a randomly selected sample from L− is greater than a random selection from E+ group.
for the groups that were found to be significantly different from the Mann-Whitney U 260 test is shown in Table 2. In this work, the CL effect size value denotes the percentage of 261 times the w f value of L− is higher than E+ when both samples are selected at random 262 from each group. It can hence be interpreted as the probability of superiority of L− 263 over E+ for a measurement. 264

265
The test-retest experiment results show that the quantitative MRI markers estimated by 266 the proposed method is repeatable. For all the markers estimated, the zero level was 267 comfortably inside the 95% confidence interval of the mean difference observed between 268 the test and retest values of the marker for the 15 ROIs (refer the Bland-Altman plots 269 in Fig. 2). Hence there are no noticeable systematic changes in the estimated markers 270 for the test and retest data [24]. considered, the medium-T 2 water pool is highly heterogeneous. It conveys information 303 on unmyelinated axons, glia, interstitial and extra-cellular matters [10].

304
Our study on MS patients has certain limitations. First, the clinical data available 305 was not of the highest quality possible due to acquisition time constraints in a clinical 306 setting. T 2 relaxometry data with a higher number of echoes and shorter echo times are 307 favorable for multi-compartment models. Second, the time gap between the first and 308 second scan was three months. A shorter interval between successive acquisitions would 309 July 4, 2018 11/14 be beneficial to study the fast evolving active lesions. 310

311
In this work we proposed a multi-compartment T 2 relaxometry model to obtain 312 quantitative estimates on brain tissues with short, medium and high T 2 relaxation times. 313 Test retest experiment results showed that the proposed method does not suffer from 314 any noticeable systematic changes in terms of the markers estimated. The study of the 315 evolution of multi-compartment T 2 relaxometry markers on 10 MS patients over a