{'fraction_non_alphanumeric': 0.0404163249144441, 'fraction_numerical': 0.020957660434991657, 'mean_word_length': 4.850099272005295, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Personalized pathology maps to quantify diffuse and focal brain damage

Background and objectives Quantitative MRI (qMRI) permits the quantification of brain changes compatible with inflammation, degeneration and repair in multiple sclerosis (MS) patients. In this study, we propose a new method to provide personalized maps of tissue alterations and longitudinal brain changes based on different qMRI metrics, which provide complementary information about brain pathology. Methods We performed baseline and two-years follow-up on (i) 13 relapsing-remitting MS patients and (ii) four healthy controls. A group consisting of up to 65 healthy controls was used to compute the reference distribution of qMRI metrics in healthy tissue. All subjects underwent 3T MRI examinations including T1, T2, T2* relaxation and Magnetization Transfer Ratio (MTR) imaging. We used a recent partial volume estimation algorithm to estimate the concentration of different brain tissue types on T1 maps; then, we computed a deviation map (z-score map) for each contrast at both time-points. Finally, we subtracted those deviation maps only for voxels showing a significant difference with healthy tissue in one of the time points, to obtain a difference map for each subject. Results and conclusion Control subjects did not show any significant z-score deviations or longitudinal z-score changes. On the other hand, MS patients showed brain regions with cross-sectional and longitudinal concomitant increase in T1, T2, T2* z-scores and decrease of MTR z-scores, suggesting brain tissue degeneration/loss. In the lesion periphery, we observed areas with cross-sectional and longitudinal decreased T1/T2 and slight decrease in T2* most likely related to iron accumulation. Moreover, we measured longitudinal decrease in T1, T2 - and to a lesser extent in T2* - as well as a concomitant increase in MTR, suggesting remyelination/repair. In summary, we have developed a method that provides whole-brain personalized maps of cross-sectional and longitudinal changes in MS patients, which are computed in patient space. These maps may open new perspectives to complement and support radiological evaluation of brain damage for a given patient.


Introduction
Neuroinflammatory and neurodegenerative diseases induce variable damage to the central nervous system, which can be quantified using magnetic resonance imaging (MRI). New advanced MRI techniques exhibit a level of sensitivity and specificity to focal and diffuse pathology that are not available through conventional MRI (Enzinger et al., 2015). Relaxometry and magnetization transfer imaging (MTI), for example, provide measures of biophysical parameters that show different sensitivity to the amount of free water, the amount of water bound to micro-and macromolecules and the amount of paramagnetic substances (i.e. iron), for review see (Helms, 2015). The combination of these MRI contrasts are sensitive and specific to the consequences of neuroinflammation and neurodegeneration, thus they are invaluable assets to perform disease diagnosis and monitoring. et al., 2015;Hilbert et al., 2018;Metere et al., 2017;Sumpf et al., 2011) and to standardize (Deoni et al., 2008;Droby et al., 2015;Helms, 2015) qMRI acquisitions, in order to improve their translation into clinical applications. Yet, how to perform individual evaluation of qMRI changes remains an open question.
In contrast to studies based on the comparison between groups, single-subject predictions (i.e. the classification of one subject as being part of one of the groups of the study, e.g. healthy vs. patient) are challenging and prone to the biases induced by the analysis method and the bias inherent to the data, which is increased in subjects with neurological pathologies (Arbabshirani et al., 2017).
In recent years, machine-learning based methods have been proposed to classify patients with pathology and healthy subjects, which provide the probability of a dataset to belong to either one of the categories. Yet, these approaches suffer from the intrinsic limitations related to (i) how feature selection is performed to differentiate the targeted classes (i.e. classical machine methods) and/or (ii) the risk of overfitting when very complex models are applied to relatively small amount of data (classical and modern machine learning methods, (Arbabshirani et al., 2017;Meyer et al., 2017;Stephan et al., 2017).
A more traditional approach to compare single-patient data to a cohort of healthy individuals is based on a voxel basis by using voxelbased morphometry (VBM) analysis (Beeson et al., 2011;Colliot et al., 2006;Kruggel et al., 2017;Muhlau et al., 2009;Sajjadi et al., 2013). Voxel-based methods rely on the non-rigid registration of the patient brain volume to a group of healthy control volumes in a common space prior to the voxel-wise comparison. The non-rigid registration can be then combined with a deformation field, which is encoded for each subject as tissue density changes in the normal space.
Despite the undeniable utility of voxel-based methods, they suffer from important limitations arising from (i) the large topological structural differences among individual brains, (ii) the necessary averaging resulting in image blurring introduced to achieve noise reduction and spatial coherence -and (iii) partial volume effects (PVE). PVE occurs when an imaging voxel contains more that one tissue type, yielding a signal equal to the weighted average of its components. The main caveat is that in regions with spatially varying intensities such as edge structures, the aforementioned effects can shift the apparent position of abnormalities detected by several voxels (Bookstein, 2001). Moreover, the interpretation of the results using single-to-group VBM analysis is problematic; especially if we consider that any statistical difference between a single subject and a group might either just reflect individual variability in neuroanatomy (instead of any alteration due to the disease under study) or designate a false positive as result of the data being sampled from a non-normally distributed population (Scarpazza et al., 2013).
In this study, we propose a new method to perform an individual analysis of qMRI maps that provides complementary information about changes in brain tissue on a single-subject basis. Different than voxelbased methods, our methodology allows a voxel-wise comparison in the subject space without the need of inter-patient registration. In short, for a given voxel in a qMRI map, we compare its scalar value with the distribution of the voxels belonging to the same tissue class and region in a group of healthy subjects and provide deviation maps from its "normative" range. To assess our method and its ability to detect subtle microstructural pathological changes within the brain, we applied our approach to a longitudinal multi-contrast dataset of relapsing-remitting multiple sclerosis (RRMS) patients, which included T1, T2, T2* relaxometry and magnetization transfer imaging.

Population
Fifteen (15) (Polman et al., 2011) and twenty (20) healthy controls (HC) (12 women/8 men, age = 31 ± 7 years) were enrolled at the Lausanne University Hospital (CHUV), Lausanne, Switzerland. A larger cohort of 65 healthy controls was used as a reference for qT1 maps, which were acquired at the CHUV in the context of another study, with the same protocol and in the same scanner applied for the current study (see Table 1).
Inclusion criteria were: RRMS diagnosis and less than five year disease duration from diagnosis; for HC: absence of previous neurological and psychiatric disease, no history of smoke and alcohol abuse and no current medication. Patients were followed up at two years (21.9 ± 2 months). At study enrollment, no patient had received corticosteroid therapy for at least three months. Immunomodulatory treatment consisted of high dosage of either interferon beta or fingolimod. Fingolimod is used as a first line therapy in Switzerland.
The ethics committee of the Lausanne University Hospital (CHUV) approved the study and all the participants gave written informed consent.
After preprocessing and qualitative evaluation of acquired MR images, 2 RRMS subjects were excluded due to motion, distortion and other artifacts seen in at least one of the contrasts used in the study (see Table 1).
The baseline and longitudinal personalized maps were computed for the 13 RRMS patients and four healthy subjects (2 women/2 men; age 34.5 ± 7.5 years), who were randomly chosen out of the control group (HC).
The healthy distribution (HD) of each contrast's parametric map used for reference was computed from the remaining 16 healthy controls (10 women/6 men; age 32 ± 10 years). An additional healthy distribution from 65 controls (42 women/23 men; age 51.68 ± 3.46 years) was used to calculate a reference for qT1 maps as stated above (HD2, Table 1). The four healthy subjects selected as a control test group for the computation of the personalized deviation map were not part of the group of controls used to calculate the healthy distribution.

Image acquisition and reconstruction
All subjects underwent MRI examinations in a 3T Magnetom Trio (Siemens Healthcare, Erlangen, Germany) using the following protocol: (1) High-resolution magnetization-prepared acquisition with gradient echo (MPRAGE); (2) 3D fluid attenuated inversion recovery (3D FLAIR); (3) 3D double inversion recovery (3D DIR); (4) Magnetization-  (Marques et al., 2010) for lesion identification and T1 relaxometry maps; (5) T2* relaxation using a prototype Gradient-Recalled-Echo (GRE) sequence with and without magnetization transfer (MT) preparation pulse (flipangle 220°, duration = 4000 ms; pulse offset = 2000 Hz; spoiler moment 25,000 μs × mT/m (Helms and Hagberg, 2009)). MTR maps were derived from the T2* data using the following formula: With M 0 and M T the images acquired without and with MT pulse, respectively. MTR was computed for each echo and averaged over all echoes. Finally, (6) T2 relaxation obtained using a prototype 2D Carr-Purcell-Meiboom-Gill (CPMG) sequence with a nonlinear inverse reconstruction algorithm that directly estimates a T2 and spin-density map from a train of undersampled spin echoes (Sumpf et al., 2011). The acquisition of T2 relaxometry maps was performed with a lower spatial resolution than the one achieved for the other MRI contrasts.
Sequence's parameters are summarized in Table 2. Visual inspection of image quality was performed in all cases.

Personalized qMRI analysis
The proposed method relies on a voxel-to-region comparison, calculating a deviation score (i.e. a z-score) between (i) the MR metric of the input image in each voxel and (ii) the MR metric distribution of the corresponding healthy tissue of the same brain region. To test the sensitivity of each parametric map to detect neuroinflammatory or neurodegenerative processes, the personalized maps were computed over the T1, T2 and T2* relaxation and MTR maps.
The condition that should be fulfilled to minimize biases is that at least one contrast should provide sufficient detailed structural information in order to quantify the concentration of each brain tissuetype, which is done to minimize PVE.
The overall methodology can be summarized as follows: (1) data preprocessing.
(2) brain segmentation into white matter (WM), gray matter (GM), cerebrospinal fluid (CSF) and a segmentation into four main lobes and deep gray matter nuclei (DGMN); (3) tissue concentration's estimation for partial volume correction and tissue's intensity characterization at each voxel; (4) segmentation and identification of lesions as WM, GM or mixed lesions; (5) computation of the data statistics (mean, variances and co-variances) of the healthy tissue intensity distribution for each tissue and region of interest; (6) computation of the deviation score for each voxel in the patient MRI based on the healthy tissue distribution previously calculated, for each qMRI contrast.
Step (4) was computed on the UNI (MP2RAGE), 3D FLAIR and DIR contrasts. An additional affine transformation to UNI image was added for the T2, T2* and MTR contrasts in order to map segmentation masks and tissue concentration maps (the masks and tissue concentration maps were computed only over the T1 maps, since they are more sensitive to differences between brain tissue types than other maps; see Section 2.3.2).
All data outputs were inspected at all processing steps to assure the accuracy of volume estimation, alignment of multi-modal MRI and fusion with lesion maps. Note that no registration is needed if only one contrast is used and the tissue characterization using the PVE is done over the same contrast. The following paragraphs provide a comprehensive description of these steps.

Data preprocessing
UNI MP2RAGE images were skull-stripped using the Morphobox prototype (see next section). T2* and MTR MR scans were corrected for gradient distortions and B1 inhomogeneities

Brain segmentation
Tissue segmentation was performed over the UNI MP2RAGE image (in patient space) using the Morphobox prototype, a software based on variational expectation-maximization tissue classification Schmitter et al., 2015). The same algorithm parcellated the brain in four main lobes: parietal, frontal, occipital and temporal (from now our regions-of-interest, ROIs), as well as the DGMN (thalamus, caudate, putamen and globus pallidus, see Supplementary Information section, SI.1, for an example of ROIs segmentation).

Tissue concentration estimation for partial volume correction
In many brain regions, especially at the border between different brain tissues, MRI intensities suffer from the so-called PVE, i.e. the imaging voxels contains more than one tissue type, yielding a signal equal to the weighted average of its components. To account for this effect in our analysis, we estimated the concentration of GM, WM and CSF in each of these neighboring voxels, assigning value computed as a mixture of tissues instead of assigning them a single tissue type. This was done by applying a PVE algorithm that was proposed by Roche and Forbes (Roche et al., 2014), which relies on continuous Markov Random Field (MRF) models, with the tissue concentration estimation formulated as a Bayesian maximum a posterior (MAP). This algorithm specifically uses a model that describes the intensity of each voxel (y) as the sum of GM, WM and CSF characteristic intensities in a given region, weighted by their respective local concentrations, with additive Gaussian noise: where μ GM, μ WM, μ CSF are the characteristic intensities of each tissue and C GM, C WM, C CSF are the concentrations of GM, WM, and CSF respectively. ε is the additive Gaussian noise (with zero mean and σ standard deviation).
In order to solve this ill-posed problem, we set three hyper-parameters to regularize the following assumptions on brain tissue intensities: (i) voxels with mild partial volume are more frequent than those with strong partial volume; (ii) tissue concentration maps are spatially smooth; and (iii) mean tissue intensities are bounded (see SI.2 The algorithm proceeds iteratively until convergence. Once the characteristic intensities μ GM , μ WM , μ CSF are initialized, the algorithm estimates the tissue concentrations C GM , C WM , C CSF (step 1) and the characteristic intensities μ GM , μ WM , μ CSF (step 2) as well as the noise standard deviation, σ. These steps are repeated until the cost function of  the algorithm converges below a defined threshold (see SI, for summary description of the PVE method). The algorithm outputs three concentration maps, one for each tissue (GM, WM and CSF) and the characteristic intensity values for GM, WM and CSF (see Fig. 1A  Since T1 maps are more sensitive to differences between brain tissue types than other maps, we performed the tissue concentration estimation on quantitative T1 maps: we initialized the mean intensities at 1350 ms for the GM and 850 ms for WM, which are the characteristic mean relaxation times according to (Marques et al., 2010). In addition, we applied the algorithm on each ROI separately to consider the variability of T1 signal over the brain (see Section 2.3.6.). We then registered the T2, T2* and MTR maps to the subjects' T1 map space (intra-subject registration) to use tissue concentrations maps estimated from the T1 contrast. To do this, we estimated an affine transformation from the T2 and T2* map space (which is the same than the MTR space) to the T1 map space using Elastix c++ library (Klein et al., 2010).

Lesion segmentation and lesion tissue identification
Brain lesions were manually counted by consensus between an expert neurologist and radiologist the in the 3DFLAIR, DIR and MP2RAGE uniform images as previously performed (Bonnier et al., 2014;Bonnier et al., 2015;Kober et al., 2012;Romascano et al., 2015). A trained technician generated manual contours for each lesion in the three different contrasts. In order to maximize the sensitivity of lesion count and volume, the masks of segmented lesions in each contrast (3D FLAIR, DIR and MP2RAGE) were then merged in a single "union lesion mask" as in (Bonnier et al., 2014;Bonnier et al., 2015;Kober et al., 2012;Romascano et al., 2015). Lesion volumes were computed and normalized by total intracranial volume as obtained using an in-house software and classified as WM, GM or mixed lesions. The lesion union mask and the ROIs masks were then registered to the T1, T2, T2*, and MTR maps.
The method used to estimate tissue concentration is robust to subtle signal intensity variations, but fails when tissue damage drastically changes MR signal like in brain lesions and their vicinity. In that case, the algorithm identifies altered WM as a mix of WM and GM (or even CSF). Thus, we a posteriori modified the concentration maps inside brain lesions, based on a single hypothesis: there is no GM inside WM. Specifically, we attributed (i) 100% WM concentration to a voxel belonging to a pure WM lesion, (ii) 100% GM concentration to a voxel belonging to a GM lesion and (iii) 50% WM and 50% GM concentration to a voxel belonging to a mixed WM-GM lesion. While in our work lesions were manually segmented, we have also integrated a software for automated lesion detection in the pipeline (Fartaria MJ, 2017a, b).
Additionally, we performed morphological operations on the white matter concentration map to apply the a-priori knowledge that there is no GM in the WM: first, we binarized the WM concentration map by thresholding it (we use a high threshold in order to keep voxels with partial volume (GM/WM border) out of the process). Then we filled the holes inside the WM binary map and finally combine the binary map with original WM concentration map (see Fig. 1C).

Computation of healthy tissue distribution for each qMRI map
T1, T2, T2* and MTR metrics vary depending on the underlying microstructure in the subjects and therefore exhibit different values  across the brain (Aubert-Broche et al., 2009;Wansapura et al., 1999). However, these parameters are considered relatively homogeneous within each cerebral lobe (Georgiades et al., 2001;Hasan et al., 2012).
For each qMRI contrast, we estimated a healthy tissue distribution for each brain lobe (frontal, parietal, occipital and temporal). We used a cohort of 65 healthy controls to compute the reference healthy tissue distribution for the T1 parametric maps and 15 controls for the other parametric maps in each brain tissue (WM, GM, and CSF) and ROI.
For each control, we first segmented the brain into WM, GM and CSF tissues as well as in four different lobes (ROIs) (Section 2.3.1) and then estimated their tissue concentrations separately (Section 2.3.2). For each control subject and ROI, only voxels with at least 95% concentration of the same tissue (WM or GM) were considered. To estimate the data's statistics (mean, standard deviation and covariances) considering intra-and inter-subjects variability, for each ROI we performed a bootstrapping analysis on the entire healthy controls cohort. In this context, we consider the data statistics (mean and standard deviation values) derived from the healthy distribution as being reasonably representative of the distribution, so that the z-score (used as a measure of deviation) shows abnormal brain tissue microstructure (see SI.3).

Deviation maps
The personalized deviation maps were computed by comparing the intensity of a single voxel in the patient scalar map with the corresponding tissue type (WM, GM or mixed) and brain ROI (frontal, parietal, temporal or occipital) value in the reference healthy distribution. The mean and standard deviation values derived from the healthy distribution of each single parameter in each ROI was considered as being reasonably representative of the distribution, so that each deviation from "healthy tissue properties" may be calculated using the zscore formula. Formally, we estimated the deviation as: where Z V is the z-score for a single voxel, Iv is the intensity of voxel v, μ is the mean of the healthy distribution over the 15 controls, and σ its standard deviation. This expression of the z-score can be extended to a more general approach taking inn consideration the different tissue concentrations. In this case, the z-core per each voxel can be re-written as: where Iv is the intensity of the voxel v in image and C GM, C WM , C CSF are the concentration of GM, WM and CSF in the voxel v, respectively. μ GM , μ WM , μ CSF and σ GM , σ WM , σ CSF are, respectively, the characteristic intensities and variances of the GM, WM and CSF tissues in the ROI where the voxel v belongs, while Cov GM/WM , Cov GM/CSF , Cov WM/CSF are respectively the covariances of the GM and WM, GM and CSF, WM and CSF in the same region (see Fig. 2)   Fig. 6. Axial view of a T1 z-score map and a T1 map in RRMS patient 1 (top row). Zoomed area showing the T1 z-score in one exemplary WM lesion (left) and the corresponding T1 map (right). The lesion appears larger in the deviation map with larger T1 z-scores in the center of the lesion and lower but still ≥3T1 z-scores in the lesion periphery/surrounding tissue.

Statistical analysis
Each deviation map quantifies the parametric variation representative of the changes occurring in an individual brain compared to a distribution of healthy controls. Thus, in order to highlight both focal and diffuse brain changes at baseline and two-years follow up, we first computed the deviation maps over the four quantitative MR contrasts (T1, T2, T2* and MTR) for both RRMS subjects and the HC group. Then, we analyzed the mean z-scores inside all voxels in all patients (i) within the lesions, (ii) in the lesions vicinity as well as (iii) in the normal appearing white matter (NAWM) and gray matter (NAGM) for all parametric maps.
The classification of "pathological" z-scores was performed based on a previous work (G. Bonnier, A. , where we MS lesions were classified into 3 groups as follows: (i) z very low (z < −2), (ii) z very high (z > 2), and (iii) z close to the HC distribution (−2 ≤ z ≤ 2), considering the continuous distribution (without distinct cluster) of lesions z-scores in each contrast: these thresholds were chosen based on the fact that in a normal distribution > 95% of the z-scores belongs to the interval [−2, 2] and that values beyond this interval reflect significant differences in patients compared to controls (P < .05).
For the longitudinal analysis, we computed difference maps (between the z-score maps obtained from the T1 map acquired at baseline and the ones obtained at 2-years follow-up) to show focal and diffuse brain changes.

Healthy tissue distribution
The assessment of the distribution of the T1, T2, T2* and MTR maps in our group of healthy controls validate the hypothesis of a reasonably normal distribution of MRI intensities in all four lobes (see Supplementary information, Section SI.3). Interestingly, the standard deviation (σ) of the intensities in the voxels within the different lobes of the healthy tissue distribution appear higher in GM than WM for all qMRI metrics (σ GM,T1 = 160 ms; σ GM,T2 = 16 ms; σ GM,T2⁎ = 22 ms; σ GM,MTR = 0.04 ms; σ WM,T1 = 60 ms; σ WM,T2 = 8 ms; σ WM,T2⁎ = 10 ms; σ WM,MTR = 0.02 ms).
Regarding the NA tissue, patients 1, 6 and 10 with a high LL but also patient 12 with a low LL (LL 12 = 0.40), showed a higher positive T1 zscore in NA white matter (NAWM) (z T1 > 0.80) than other patients. All patients showed a low z-scores in NAGM, with patients 1 and 13 have the highest deviation (z T1 > 0.40) and patients 3 and 4 have the lowest deviation (z T1 < 0.10). The highest z-scores were observed around and between lesions (Fig. 5). Also, MS lesions in the T1 deviation maps often appeared to be larger than the manually segmented lesions (Fig. 6). This was most probably due to areas of increased T1 relaxation times surrounding MS lesions that were not perceived as lesion tissue by our neurologist and radiologist experts.
Finally, the z-score map inside some large lesions revealed presence of "core" of increased T1 relaxation time where the z-score is even higher than in the rest of the lesion (Fig. 7).

Longitudinal analysis
In patients, some subcortical areas, predominant in the fronto-parietal lobes, showed a concomitant increase in T1, T2, T2* z-scores and Fig. 10. Axial difference between baseline (time point 1) and 2 years follow up (time point 2) of T1, T2, T2* and MTR z-score maps in patient 6 (top), patient 1 (middle) and control 2 (bottom). The difference maps give an overview of the changes occurring in the brain tissue of patients 6 and 1, and control 2. decrease of MTR z-scores over two years ( Fig. 10A and C). In the lesion periphery, some areas e exhibited a decrease in T1/T2 and slight decrease of T2* z-scores ( Fig. 10A and B). Besides, other areas showed a decrease in T1, T2 (T2*) and a concomitant increase in MTR (Fig. 10B).

Discussion
We have developed a novel approach to provide a personalized quantification of brain damage using a voxel-wise comparison between quantitative parameters of tissue integrity in a patient and the distribution of the same parameters in the same region and tissue in a population of healthy controls.
To test the sensitivity and robustness of our method in the detection of subtle pathological brain changes, we have applied it to a longitudinal dataset of early-stage RRMS patients, who were imaged using four quantitative and semi-quantitative contrasts (T1, T2, T2*, MTR). Those contrasts/metrics were chosen because their combination has been previously proven to increase sensitivity and specificity to pathological brain changes related to inflammation and degeneration compared to single contrast/metric approaches (Bonnier et al., 2014;Bonnier et al., 2015;Granziera et al., 2015;Romascano et al., 2015).
When both RRMS and control subjects were analyzed at baseline, the proposed method was able to detect deviation from "normative ranges" (computed as healthy tissue characteristic intensity obtained from healthy subjects) for RRMS subjects, while no significant deviations were found for our control test group. These differences were especially remarkable within MS lesions and their vicinity but were also sometimes evident in the remaining normal-appearing tissue.
As expected from previous group studies (Bonnier et al., 2014;Bonnier et al., 2015;Deoni et al., 2008;Neema et al., 2007), T1 relaxation deviation maps computed in RRMS patients showed a global increase in T1 relaxation time compared to the healthy control population, which suggests the loss of tissue integrity (i.e. axonal, myelin and cell damage) and/or extracellular water accumulation.
In MS lesions, T1 z-scores were in general increased and some lesions showed the peculiar microstructure reported in postmortem studies i.e. a gradient of damage from a central core with more pronounced tissue loss (i.e. more prolonged T1 relaxation times compared to healthy tissue) and areas surrounding it with less prolonged T1 relaxation times (Figs. 6, 7) (Bruck et al., 2002;Lassmann et al., 2001;Lucchinetti and Bruck, 2004;Metz et al., 2014).
In the normal appearing tissue, the highest T1 z-scores were measured in the regions surrounding MS lesions, as previously reported (Bonnier et al., 2014;Davies et al., 2007;Parry et al., 2002;Seewann et al., 2009;Vrenken et al., 2010;West et al., 2014), (Figs. 3, 5). Increased T1 relaxation times point at a loss of tissue microstructural damage and/or to an increase in iron content in the tissue (Bonnier et al., 2017;Bonnier et al., 2014;Bonnier et al., 2015;, which are visible also in the absence of macroscopic global or localized changes in brain volume (Bonnier et al., 2017;Bonnier et al., 2014).
The remaining normal appearing tissue showed mild T1 increases compared to the healthy population of references and interestingly, NAGM showed lower deviation in T1 relaxation times compared to healthy tissue than NAWM. This may be due to the presence of a lower degree of tissue damage in GM than WM at early disease stages, which is consistent with previous studies (Davies et al., 2007;Klaver et al., 2013;Muhlau et al., 2013). Yet, this may also be due to higher PVE, which induces a higher variability when computing the healthy GM distribution compared to the WM.
As well, and as expected from previous literature (Parry et al., 2002(Parry et al., , 2003Rovira et al., 2013), the evaluation of T1 z-scores deviations in MS patients showed a diffuse microstructural pathology that is not related to lesion load. Therefore, the personalized maps of tissue damage -obtained with the proposed methodology-may help to better depict the global damage provoked by MS disease, beyond the areas of focal lesions.
On the other hand, the deviation maps obtained from T2, T2* and MTR maps showed brain tissue alterations, which were located in the same regions as the ones obtained from T1 maps (Fig. 8). Yet, while the T2 and T2* deviation maps highlighted more diffuse changes in patient white matter and gray matter tissue, the MTR deviation maps mainly revealed focal changes in lesions (Fig. 8). Lesion and normal-appearing tissue z-score deviations were positive for T2, T2* maps and negative for MTR maps (Fig. 9), which is compatible with (i) tissue degeneration/loss (in the presence of concomitant and strong increase in T1) or (ii) inflammatory phenomena (in the presence of less pronounced increase in T1 and/or concomitant increase in T2/T2*).
Nonetheless, the z-scores measured in T2, T2* and MTR maps were less important than the ones revealed by T1 maps: this may be due to the lower spatial resolution of the T2, T2* and MTR maps, which translates in an increase of partial volume effect, but also to a higher variability of the healthy distribution compared to the pathological effect on the T2 and T2* relaxation times.
Longitudinal analysis of parametric deviation maps showed no changes in control subjects but areas of z-score alterations in patients. Specifically, in RRMS patients, we have evidenced regions with a concomitant increase in T1, T2, T2* z-scores and decrease of MTR z-scores, a behavior suggesting brain tissue loss over time (Bonnier et al., 2017;Bonnier et al., 2014;Deoni, 2011;Filippi and Agosta, 2007). Also, in the lesion periphery, we could observe areas with decreased T1/T2 and slight decrease of T2* indicating iron accumulation (Dusek et al., 2013). Moreover, in other areas, we assessed longitudinal decrease in T1, T2 (T2*) and a concomitant increase in MTR occurs, suggesting remyelination/repair (Mallik et al., 2014), Fig. 10.
In summary, we have developed and assessed a method that provides a whole-brain personalized information of brain tissue damage and showed its sensitivity to diffuse and focal changes in MS patients, both cross-sectionally and longitudinally. One advantage of the proposed method compared to conventional voxel-based morphometry is that the maps are computed in patient space, therefore bypassing all limitations of inter-subject registration (Klein et al., 2010). Another advantage is the fact that it may be easily extended to other type of images, provided that a large number of healthy subjects are also imaged with the same settings. Also, the proposed approach may be easily extended to other neurological diseases since the applied metrics are sensitive to neurodegenerative and neuroinflammatory changes also in other neurological disorders such as mild cognitive impairment , neuroHIV infection (Granziera et al., 2013) and migraine .
The proposed method may still benefit of a number of targeted improvements. Indeed, the healthy cohort of reference presently consists of 65 and 15 subjects, who are reasonably but not perfectly representative of the healthy distribution of each metric. Because of this, the presented results require confirmation using a larger reference cohort of healthy controls, which may reduce the variability of the intensity distribution, increase the statistical power and improve sensitivity of the comparison.
In addition, the T2 mapping method that we have applied in this work models a mono-exponential T2 signal decay (Sumpft et al., 2011), which may be violated by the presence of stimulated echoes. To minimize this issue, we have discarded the first acquired echo (at TE = 9 ms). Despite that, however, the signal sampled with the remaining echoes (in the range from 20 to 170 ms) may be still subjected to model violations by the different T2 contributions in individual brain voxels (MacKay, Whittal, et al. 1994, MacKay, Laule, et al. 2006, Barta, et al. 2015. Indeed, for a precise measurement and quantification of the different T2 components of the brain tissue, a much longer acquisition and a more complex data modeling are necessary (MacKay, Whittal, et al. 1994, MacKay, Laule, et al. 2006, Barta, et al. 2015. Nevertheless, the applied T2 mapping measurement has shown high sensitivity to brain pathological changes in neuroinflammatory and neurodegenerative diseases (Bonnier et al., 2017;Bonnier et al., 2014;Bonnier et al., 2015; as well as high precision and reproducibility (T2 variation of < 0.5 ms were reported in scan-rescan measurements and < 4 ms across subjects) (Hilbert, et al. 2018) and relatively short scan times, all characteristics that render it suitable for clinical studies. Furthermore, the approach we propose uses the information provided by different contrasts and metrics (i.e. T1, T2, T2* and MTR) to identify pathological changes in brain tissue microstructure and this is complementing -to a certain extent -the missing information derived from a mono-exponential fitting of the T2 decay.
Lastly, future studies should also further explore the optimal cut-off or range to detect neurodegenerative and inflammatory processes.
In conclusion, we have provided a method to calculate personalized deviation maps in MRI, which may open new perspectives to complement and support clinical decisions. Future work will focus on increasing the healthy control population, integrating automatic lesion detection tools in the pipeline and in the identification of the optimal threshold for pathological changes in each contrast.