Improved myelin water imaging using B1+ correction and data-driven global feature extraction: Application on people with MS

Abstract The predominant technique for quantifying myelin content in the white matter is multicompartment analysis of MRI’s T2 relaxation times (mcT2 analysis). The process of resolving the T2 spectrum at each voxel, however, is highly ill-posed and remarkably susceptible to noise and to inhomogeneities of the transmit field (B1+). To address these challenges, we employ a preprocessing stage wherein a spatially global data-driven analysis of the tissue is performed to identify a set of mcT2 configurations (motifs) that best describe the tissue under investigation, followed by using this basis set to analyze the signal in each voxel. This procedure is complemented by a new algorithm for correcting B1+ inhomogeneities, lending the overall fitting process with improved robustness and reproducibility. Successful validations are presented using numerical and physical phantoms vs. ground truth, showcasing superior fitting accuracy and precision compared with conventional (non-data-driven) fitting. In vivo application of the technique is presented on 26 healthy subjects and 29 people living with multiple sclerosis (MS), revealing substantial reduction in myelin content within normal-appearing white matter regions of people with MS (i.e., outside obvious lesions), and confirming the potential of data-driven myelin values as a radiological biomarker for MS.

rate of water trapped between myelin sheaths and the relaxation rates of the intra-/extracellular water pools.Prominent techniques within this category include gradient and spin echo (GRASE) ( Prasloski, Rauscher, et al., 2012), multicomponent driven equilibrium single-pulse observation of T 1 and T 2 (mcDESPOT) ( Deoni et al., 2013;Zhang et al., 2015), T 2 magnetization prepared imaging ( Nguyen et al., 2012;Oh et al., 2006), and magnetic resonance fingerprinting ( Y. Chen et al., 2019;Nagtegaal, Koken, Amthor, & Doneva, 2020).To fully utilize the differences in T 2 relaxation times and avoid biases caused by field inhomogeneities, approaches based on spin echo protocols are favorable.Signal processing then involves multicomponent T 2 (mcT 2 ) analysis ( Does, 2018;Whittall et al., 1997), where one maps the relative fraction of the fast-relaxing T 2 component (<40 ms) corresponding to water bound between myelin sheaths, and the longer T 2 component (40-200 ms) corresponding to intra-/extracellular water ( MacKay & Laule, 2016).The myelin water fraction (MWF) serves as an indirect measure of myelin content and is calculated as the ratio of the T 2 spectral content corresponding to myelin water to the total area under the T 2 spectrum.
The most efficient protocol for mcT 2 mapping in vivo is 2D (i.e., multislice) multiecho spin echo (MESE), offering high-sensitivity T 2 relaxation times at clinically relevant scan times of 5-8 minutes.MESE signals, however, are contaminated by stimulated and indirect echoes ( Hennig, 1988), thereby deviating from ideal multiexponential behavior.Various strategies have been proposed to mitigate these effects, e.g., shifting to 3D MESE acquisitions ( Dvorak et al., 2020;Meyers et al., 2009;Prasloski, Mädler, et al., 2012), adjusting the refocusing slice width to be significantly larger than the excitation slice profile ( Faizy et al., 2016;Kumar et al., 2016;Pell et al., 2006), or using complex crusher gradient schemes ( Kolind et al., 2009;Mackay et al., 1994).A more effective solution is the utilization of advanced algorithms such as the extended phase graph (EPG) method, which iteratively traces the evolution of multiple spin populations throughout an MESE echo train by incorporating the stimulated echoes into their signal model ( Hennig, 1988;Lebel & Wilman, 2010).Another efficient approach is the echo modulation curve (EMC) algorithm, which employs Bloch simulations to comprehensively account for all coherence pathways arising during the echo train, thereby faithfully reproducing both stimulated and indirect echoes ( Ben-Eliezer et al., 2015, 2016;McPhee & Wilman, 2017;Radunsky et al., 2021).
Once a model is chosen for resolving the bias due to stimulated echoes, the fundamental approach to multicomponent analysis involves the inversion of the signal acquisition process.This inversion aims to identify the signals originating from distinct cellular water pools and extract the relative signal intensities associated with each of these pools.This inverse problem is inherently ill-posed and poses a significant challenge, even when incorporating various regularization techniques.Specifically, this inversion process exhibits nonuniqueness in the solutions space and a high susceptibility to noise ( Graham et al., 1996;Whittall & MacKay, 1989).Consequently, despite its paramount relevance of myelin mapping for clinical applications, the field of MWI still lacks a gold standard which can effectively address these complexities.
In this study, we utilize a new paradigm for mcT 2 fitting, which effectively mitigates issues related to ambiguity to produce reliable MWF maps.The novel aspect of this technique is a preprocessing step which performs statistical analysis of the signals from the entire WM in order to identify a set of multicomponent configurations (termed mcT 2 'motifs') which best describe the examined tissue.These motifs are then used as basis functions for a regularized non-negative least square (RNNLS) mcT 2 fitting of the signal within each voxel ( Omer et al., 2022).To further enhance its accuracy, our approach relies on the EMC signal model, which incorporates the exact pulse-sequence scheme and scan parameters to address the existence of stimulated and indirect echoes.This integration ensures the provision of accurate and reproducible T 2 values that remain consistent across scanners and scan settings ( Ben-Eliezer et al., 2015, 2016;Radunsky et al., 2021).Important methodological improvements are introduced in this study including accounting for transmit field (B 1 + ) inhomogeneities, the use of entropy-based regularization constraint, and enforcing pseudo-orthogonality among mcT 2 motifs.Validations are presented on a numerical phantom at varying SNRs, and on a physical three-compartment phantom having a unique design which offers ground truth values.Performance of the data-driven approach is also evaluated vis-à-vis conventional RNNLS processing.Clinical applicability of the new approach is demonstrated on healthy subjects and people with MS.

Conventional EMC-based mcT 2 analysis
For readers' convenience, we provide a concise overview of the conventional, nondata-driven approach to mcT 2 fitting of MESE signals using the EMC signal model and RNNLS fitting.A comprehensive description of this method is also available in Omer et al. (2022).
The MESE signal from each voxel represents a superposition of signals from several subvoxel (cellular) water pools.Due to contamination by stimulated and indirect echoes, attempting to model this signal using a multiexponential decay approach will lead to inaccuracies.To address this issue, we employ the EMC algorithm, which generates realistic T 2 decay curves by accurately simulating the exact pulse-sequence scheme, radiofrequency (RF) and gradient pulse shapes, and timing diagram of the 2D MESE protocol, used for image acquisition.The signal from each voxel can then be modeled as: where s ∈R ETL is the experimentally acquired signal, is the relative fraction of each subvoxel water pool, ETL is the echo train length, and N T 2 is the number of cellular components.Extraction of w is typically done by solving an RNNLS minimization problem of the form: where λ Tikh , λ L 1 ≥ 0 are Tikhonov and L 1 regularization terms.Tikhonov regularization is added to favor smoother solutions, while the L 1 regularization is added to promote sparse solutions.Due to the large space of possible mcT 2 combinations, more than one solution can match each experimental signal.Adding to that the ambiguity caused by noise, solutions of this system of equations tend to be highly unstable and, in many cases, converge toward a wrong local minimum.

Data-driven EMC-based mcT 2 analysis
In order to address the inherent ill-posed nature of the inverse problem, we propose the implementation of a data-driven preprocessing stage aimed at identifying distinctive mcT 2 motifs that best encapsulate the characteristics of the WM tissue.This preprocessing step is carried out prior to the optimization process expressed in Eq. ( 2).The underlying principle of this approach is based on one key assumption: that within the WM, there exists a finite set of microstructural configurations, each of which corresponds to a specific mcT 2 spectrum.A flowchart of the algorithm is presented in Figure 1.
The input for the data-driven algorithm consists of MESE data, the exact pulse sequence scheme and acquisition parameters, and a WM mask.The algorithm consists of eight sequential steps.
Step #1 involves the creation of a single-T 2 EMC dictionary of signals, spanning N T 2 logarithmically spaced T 2 values ranging from 10 to 800 ms.In step #2, a theoretical dictionary of all possible mcT 2 motifs is constructed by combining series of N Comp single-T 2 signals, each with a relative fraction f n ∈ 0..1 [ ] using a fraction resolution of Δf = 1/N frac , where N frac is the number of discrete values that f n can assume.Each simulated mcT 2 dictionary element ( is expressed as where N comp is the number of compartments.Importantly, the number of compartments in the theoretical dictionary does not limit the number of compartments in the final mcT 2 spectrum for each voxel, which will comprise a linear combination of mcT 2 motifs.An important substage of step #2 involves pruning nonphysiological configurations from the mcT 2 dictionary, including mcT 2 motifs that lack a short T 2 component (<40 ms); motifs whose short T 2 component fraction exceeds 30%; and motifs that correspond to a single-T 2 value that is outside the range of T 2 values that exist in the WM region.This single-T 2 value of an mcT 2 motif is calculated by fitting the motif's signal decay curve to a dictionary of single-T 2 signals.Stage #3 consists of correcting for transmit field (B 1 + ) inhomogeneities, and is elaborated in Section 2.3.In steps #4-5, statistical correlation is performed between each dictionary element ( i and every voxel signal s j , by calculating a "cost" α i, j , defined as the L 2 norm difference between each pair.A regularization term β is added to this cost value to penalize motifs with high entropy, retrieving the simplest spectral features and reducing overfitting ( L. Chen et al., 2002;Widjaja & Garland, 2005): where f i,n is defined according to Eq. (3).Prior to selecting the final basis set of motifs, we define a similarity criterion between a dictionary element (( i ) and voxel signal (s j ).This is used to apply cost clipping ( Taal et al., 2011) for the purpose of reducing potential bias from motifs which are significantly different from the experimental data.The similarity threshold is denoted as ξ MV and defined as where δ MV is set empirically.Each ( i and s j pair is then considered similar if it obeys: In this process, α i, j is constrained by α i, j → min(α i, j , 5ξ MV ), and κ i, j is updated accordingly.In step #6, the κ i, j scales are normalized to [0…1] range according to In step #7, the scores κ i, j of each dictionary element ( i are summed across all voxels to yield a global score Κ i that expresses how well each motif matches the examined tissue: (10) Finally, in step #8, a finite subset of motifs is selected, which will be used to separate each voxel's signal into its underlying components.The selection process is based on two criteria: the highest Κ i score and a pseudoorthogonality constraint which maximizes the differences between motifs in the final subset.A full step-by-step description of the selection process is detailed in Algorithm 1.

Algorithm 1
Input: Set of all mcT 2 motifs sorted according to their score Set of all the voxels in the analyzed ROI Group of voxels that are similar to motif i (see Eq. ( 8)) Output: Algorithm: Initialize D to be an empty set.
Add the mcT 2 motif with the highest score to the set of selected motifs.
Remove from the group of all voxels V !, voxels that are similar to motif ( 1 4: for i = 2, 3…,N motifs : Loop over all remaining mcT 2 motifs, sorted according to their score Check whether motif i does not have the same single-T 2 value as any of the other motifs already in D Check whether this motif holds new information, i.e., is it similar to voxels that were not yet accounted for by motifs already in D 7: D → D ∪ ( i Add the motif to the set of selected motifs 8: Remove from the group of all voxels V !, voxels that are similar to motif i 9: end if 10: If all voxels have been accounted for (i.e., V ! is an empty set) 12: break for loop Break for loop and terminate the algorithm 13: end if 14: end for After selecting the final set of motifs D, the signal is modeled as a linear combination of mcT 2 signals: where s ∈R ETL is the experimental signal, W ∈R D is the unknown vector of weights of the elements in D, and |D| is the number of elements in D. The mcT 2 fitting problem now converts to solving for the unknown vector W through a standard RNNLS optimization procedure: which is similar to Eq. ( 2), albeit with a modified encoding operator that has been learned from the tissue being analyzed.Notably, the result obtained from Eq. ( 12) is W, although our objective is to derive vector w.Recalling that each motif is defined by a weights vector f (Eq.( 3)), and denoting the matrix of all weights as F∈R N T 2 × D , the final T 2 spectrum at each voxel will be Finally, MWF values are calculated from the T 2 spectrum of each voxel as the relative energy between 0 and 40 ms, and the energy of the entire spectrum.

Correcting for transmit (B 1 + ) field inhomogeneities
The extensive utilization of RF refocusing pulses in MESE signals can introduce bias due to transmit field inhomogeneities.To address this, the data-driven algorithm has been added with a preprocessing stage designed to estimate the B 1 + profile and subsequently correct the experimental data.As each theoretical EMC decay curve d i depends on both T 2 and B 1 + , the same holds for each mcT 2 motif which can now be expressed as ((B 1 + ).Thus, each motif ( transitions from a 1D to 2D by adding a B 1 + dimension, discretized across N B 1 + values in the range 80-120% (where 100% represents a fully homogeneous field).
The B 1 + correction procedure involves three stages.
First, the initial B 1 + profile is calculated for each voxel j by finding the dictionary motif ((B 1 + ) that has the lowest L 2norm difference to the experimental signal s j .This is done using an exhaustive search over all dictionary ele-ments and produces an initial solution for the transmit field map B 1, n= 0 + .In the second stage, the field map undergoes iterative spatial smoothing within a region surrounding each voxel by minimizing the following cost function: Here, B 1,n + and B 1,n+1 + denote the transmit field profile at iteration n and n + 1, respectively, and B 1 + denotes the simulated range of N B 1 + transmit field profile values.The cost function employed in this optimization consists of two components: first, a likelihood term that is responsible for finding the B 1 + value that best aligns with the data from voxel s j ; and second, a prior that imposes spatial smoothness within a 2D kernel using L 1 norm.The utilization of the L 1 norm in this process has been shown to enhance resilience against noise and outliers across various signal processing applications ( Brooks et al., 2013;Markopoulos et al., 2014).μ represents the regularization weight, N k denotes all voxels within the 2D kernel surrounding voxel s j , and N k is the number of voxels in N k .The iterative process is terminated either when there is no change in the B 1 + value between iterations or when the number of iterations exceeds the predefined limit of N iter = 200.The resulting transmit field map is then denoted as B 1, opt + .Finally, the B 1, opt + map is used to correct signal s j in each voxel.Taking the mcT 2 motif with the optimal B 1 + value ((B 1,opt + ) , the corrected signal s j,cor is calculated as where ( denotes the "homogeneous" mcT 2 motif, corresponding to B 1 + value of 100%.

Validation on a numerical phantom
To evaluate the suggested method, numeric simulations of a 2D MESE protocol were performed on a Shepp-Logan phantom, using matrix size = 90 x 90, ETL = 11, echo time (TE) = 12 ms, interecho spacing = 12 ms, and bandwidth = 200 Hz/Px.Simulations were repeated for five tissue types, varying by the number of compartments (N Comp = 1,2,3), the myelin fractions, the relaxation times, and relative fractions of the intra-/extracellular water pools, and the B 1 + field profile.Detailed description of simulated mcT 2 configurations is delineated in Figure S1.To benchmark the myelin mapping algorithm across different SNRs, Rician noise was introduced to the simulated signals at SNRs of 500, 300, 200, 100, and 50, defined as the ratio between the first echo amplitude and the standard deviation of the noise.mcT 2 fitting was performed according to the algorithm described in Figure 1, starting with generating an mcT 2 dictionary.The number of possible configurations (N mcT 2 ) grows exponentially with N T 2 , N frac , and N B 1 + , which are used to construct the dictionary, and can be calculated using combinatorics.For example, for a choice of two compartments, this number will be which comes up to 3,404,700 elements for N T 2 = 200, Δf = 0.05, N B 1 + = 9 (B 1 + = 80 : 5 : 120%).The regularization weights λ Tikh , λ L 1 , and λ ent were determined through an exhaustive search in the range of 0-10, optimized for maximal MWF accuracy.The importance of choosing the right λ Ent is further demonstrated in Figures S2 and S3, for SNRs of 500 and 100, respectively.The minimization problem in Eq. ( 12) was solved using MATLAB's (Mathworks, version 2022b) Quadratic programming (see Appendix A for detailed description).To assess the algorithm's stability, mean absolute error was computed as a function of λ Tikh and λ L 1 values, for the five tissue types in the numeric phantom.Optimal reconstruction parameters for the conventional and datadriven approaches, i.e., those which yielded the highest accuracy, are detailed in the 3rd and 6th columns of Table 1, respectively.

Validation on a physical phantom
An mcT 2 phantom was prepared using MnCl 2 solutions at concentrations of 0.11, 0.15, and 0.6 mM, producing T 2 relaxation times of 80, 60, and 20 ms, respectively.Nine tubes, each with a volume of 70 ml, were used as containers of a physical multicompartment phantom, featuring three distinct internal compartments.The first compartment consisted of background solution with T 2 of 80 ms.Varying number of 3 and 5 mm capillary tubes were then inserted into the container tubes, filled with 60 and 20 ms solutions, and serving as the 2nd and 3rd compartments.This resulted in nine relative fractions of the short T 2 (20 ms) compartment, equal to 0%, 0%, 3.8%, 7.1%, 11.1%, 14.3%, 18.5%, 21.7%, and 26.2%.Full description of the nine different multicompartment configurations is delineated in Table S1.
A key aspect of these scans is the huge voxel size = 25 x 25 mm 2 , ensuring that each tube was fully encapsulated within a single voxel.This unique setup enabled the generation of genuine, experimental mcT 2 signals, with known ground truth fractions of the short T 2 component.To generate a large number of mcT 2 signals, each low-resolution scan was repeated four times with varying slice offsets.Lastly, to assess the interscan repeatability of the fitting process, two of the tubes having fractions of 7.1% and 18.5% were scanned twice, resulting in an overall number of 440 separate signals (10 slices x 4 offsets x 11 tubes).
Postprocessing and statistical analysis of the phantom data included the determination of ground truth MWF values based on the relative area of the short T 2 tubes, derived directly from the phantom geometry.The data-driven analysis utilized data from all tubes, creating a diverse dataset, with a range of single-T 2 values that is similar to what is expected in vivo.Mean and standard deviation (SD) of the MWF were calculated for each tube, and Pearson correlation and root mean square difference (RMSD) were computed to assess the agreement with the ground truth MWF values.Interscan repeatability was assessed by calculating RMSD and correlation coefficient, and performing Bland-Altman analysis for the two tubes which were scanned twice.Detailed list of all postprocessing parameters which gave the highest accuracy for both approaches is outlined in Table 1, 4th and 7th columns.

MWF mapping of healthy subjects
A group of 26 healthy subjects (39.2 ± 5.5 y/o, 15 males) was scanned on a 3T Prisma MRI scanner (Siemens Healthineers), after signing an informed consent and under Helsinki approval by Sheba Medical Center Prior to MWF fitting, MESE data were denoised to enhance SNR using the MP-PCA algorithm ( Stern et al., 2022), following a previous study that demonstrated the benefits of denoising for multicomponent relaxometry fitting ( Does et al., 2019).Extraction of WM mask was performed on the MPRAGE images, followed by registration to MESE image space using Freesurfer tools ( Fischl, 2012;Reuter et al., 2010).MWF maps were generated using the conventional and data-driven techniques, using the postprocessing parameters detailed in Table 1, 5th and 8th columns, and a fixed kernel size of 15 x 15 mm 2 for B 1 + correction.Single-T 2 maps were also calculated based on the mcT 2 spectrum within each voxel by (i) generating a single-T 2 signal using a linear combination of all T 2 components (see Eq. ( 11)), and (ii) performing single-T 2 fitting using the EMC algorithm ( Radunsky et al., 2021).Reproducibility of the data-driven MWF fitting algorithm was tested using scan-rescan available data from 22 out of the 26 healthy subjects.The time difference between the two scans was 30 ± 13 days. .MWF maps were generated using the same parameters as for the healthy subjects.T 2 maps were calculated from the mcT 2 spectrum at each voxel, following a procedure described above for the healthy subjects.

Statistical analysis
Mean and SD of MWF values were calculated using the data-driven approach for six manually segmented 2D normal-appearing WM (NAWM) regions.These included the genu of the corpus callosum (GCC), splenium of the corpus callosum (SCC), frontal lobe, temporal lobe, occipital lobe, and the entire WM (segmented automatically using Freesurfer software).Illustration of the segmented ROIs is shown in Figure S4.Scan-rescan measurements were assessed for repeatability by calculating correlation coefficient and performing Bland-Altman analysis.
To compare MWF values between healthy subjects and people with MS, a two tailed unpaired t-test with significance level of α significance = 0.001 was performed.Bonferroni multiple comparisons correction was applied by dividing α significance by the number of tests, which was 6 in this case.Classification of people with MS vs. healthy controls was performed based on MWF values at each normal-appearing region, and receiver operating characteristic (ROC) curves were generated for each ROI, followed by calculating the area under the curve (AUC) as a metric of classification performance.

RESULTS
The performance of the data-driven fitting of numerical phantom data is shown in Figure 2 along with a comparison with the conventional RNNLS approach.Mean absolute errors of 0.2%, 0.5%, 0.7%, 1.2%, and 1.8% were produced by the data-driven approach for SNRs of 500, 300, 200, 100, and 50, respectively.Conventional fitting, in contrast, showed consistently higher mean absolute errors of 2.4%, 2.5%, 2.5%, 2.8%, and 3.7%  across the tested SNR values.B 1 + bias field correction using the data-driven approach also showed high correlation to ground truth values with mean absolute errors of 0.1%, 0.1%, and 0.4% for SNRs > 100, while the conventional approach struggled to estimate B 1 + accurately and produced mean absolute errors of 6.0%, 5.1%, and 4.7%.At SNRs ≤ 100, higher deviations were observed for the B 1 + maps using both approaches with mean abso-lute errors of 2.8% and 5.5% for the data-driven technique and 5.3% and 6.1% for the conventional fitting at SNRs 100 and 50, respectively.In Figure 3, stability maps are provided, depicting conventional and data-driven fitting performance across λ Tikh and λ L 1 regularization values and different SNRs.The datadriven approach exhibits consistently small mean absolute errors over a wide range of L 1 and Tikhonov regularization weights for all tested SNRs.In contrast, the conventional approach demonstrates large mean absolute errors throughout the entire range of regularization settings.
The design of the three-compartment physical phantom is shown in Figure 4A-B.Fitted MWF versus ground truth fractions using the data-driven technique are shown in Figure 4C, exhibiting high linear correlation (r) and small absolute error of 0.83 ± 0.51%.Similar analysis of conventional mcT 2 fitting is shown in Figure S5.
Scan-rescan analysis of the physical phantom is presented in Figure 5A-B for two tubes with MWF fractions of 7% and 18.5%.The data-driven approach    exhibits small RMSD of 0.76% and 0.22% for the 7% and 18.5% tubes, respectively.Further statistical analysis is presented using the Bland-Altman plots in Figure 5C-D, demonstrating a good agreement between the scans with difference of 0.40 ± 0.66% and 0.15 ± 0.16% for the 7% and 18.5% tubes, respectively.Similar analysis of conventional mcT 2 fitting is shown in Figure S6.
T 2 -weighted images, T 2 maps, and MWF maps from three representative healthy subjects are presented in Figure 6.No correlation was observed between the MWF and T 2 values, i.e., higher MWF values do not necessarily indicate lower T 2 values as can be seen in the GCC and SCC regions.In vivo repeatability of MWF maps was assessed using scan-rescan data from 22 healthy subjects, summarized in the Bland-Altman and correlation plots in Figure7.The data-driven approach yielded realistic MWF values and demonstrated good repeatability, with a correlation coefficient of 0.91 and RMSD of 1.08%.Similar analysis of conventional mcT 2 fitting for healthy subjects is shown in Figures S7-S8.
Figure 8 shows FLAIR images, T 2 maps, and MWF maps of three representative people with MS.The advantage of the data-driven approach is clearly demonstrated in this figure, showing how quantitative maps reveal subtle changes which are indicative of inflammation and demyelination within NAWM, and which are not visible in the qualitative FLAIR images.Similar quantitative maps produced using the conventional mcT 2 fitting are shown in Figure S9.
MWF values for specific brain ROIs were extracted for both healthy subjects and people with MS, and then used to classify subjects between these groups.Box plots of the extracted values are illustrated in Figure 9, along with classification ROC curves.The data-driven approach produced a significant difference (p-value < 0.0001) in mean MWF between healthy subjects and people with MS across all ROIs, with a relative reduction in MWF ranging from 20% to 38%.Full numeric values per ROI are given in Table S2.The data-driven approach yielded consistently higher AUC for all tested regions.Similar analysis using the conventional mcT 2 fitting is shown in Figure S10 and Table S3.Imaging Neuroscience, Volume 2, 2024

DISCUSSION
The significant ambiguity within mcT 2 search space poses a substantial obstacle for achieving reliable MWF values.This study introduces several improvements to a recently developed MWI technique based on mcT 2 analysis of MESE data, with corrections for transmit field inhomogeneities.The method begins by identifying global mcT 2 motifs of the entire tissue, which are then used for localized signal analysis at each voxel, while accounting for local variations in the B 1 + field.As an initial step, a comprehensive dictionary of mcT 2 signals is generated, tailored to match the precise pulse-sequence parameters of the MESE protocol.Extraction of specific pulse-sequence parameters is necessary in order to use the EMC model of MESE signals, which stands at the basis of the data-driven MWF mapping algorithm.On Siemens scanners, this is done via the IDEA software.Accessing these parameters requires proprietary scanner code and a research agreement with the corresponding vendor.
The identification of tissue-specific mcT 2 motifs significantly reduces the number of potential solutions, thereby alleviating the intrinsic ill-posed nature of mcT 2 analysis.The primary objective of this stage is to keep the most relevant motifs-either having good correlation to a large number of voxels or very good correlation with a limited number of voxels, like lesions.The ensuing set of motifs constitute a pseudo-orthogonal basis set with maximum information about the tissue.This strategy mitigates the risk of converging toward less physiological solutions, even if they might demonstrate higher fitting accuracy at the voxel level.Although this paper primarily focused on 2D MESE sequences which can be run at clinical scan times of 5-8 minutes (depending on image resolution), the data-driven approach can also be adapted for use with 3D MESE sequences using dictionaries of exponentially decaying signals.Recently, another datadriven approach has been suggested ( Piredda et al., 2022), where different regression models are trained using single relaxometry measurements (such as single T 1 and single T 2 ) to predict MWF values.While this approach focuses on establishing a generalized model for MWF prediction across acquired data, the method presented herein is based on learning specific characteristics of the WM of each individual.Similarly, the proposed approach has the potential to extend its applicability by conducting the data-driven analysis on a collective group of individuals.However, such generalization should be approached cautiously, considering various factors such as age, gender, and other cofactors that may affect myelination patterns.
The number of compartments used to generate the simulated mcT 2 dictionary may be either two, three, or even higher if sufficient computational resources are available.Here, we conducted tests with both two-and three-compartment dictionaries and determined that a two-compartment dictionary is sufficient and does not introduce higher errors in the measured MWF, while allowing to perform analysis on a standard PC.Importantly, the number of compartments in the mcT 2 dictionary does not impose any constraints on the final number of components in the T 2 spectrum, as there is no inherent limit on the number of mcT 2 motifs used to represent the signal in each voxel (see Eq. ( 11)).
The range of single-T 2 values used to filter nonphysiological mcT 2 motifs can exhibit variations among subjects and may be wider in tissues containing pathologies such as MS lesions.It is, therefore, crucial that the range and dynamic resolution of T 2 values in D be sufficiently dense within the physiological range of T 2 values to construct an mcT 2 dictionary D that will faithfully characterize the tissue.Various choices of N T 2 have been reported in the literature ( Kumar et al., 2018;Nagtegaal, Koken, Amthor, De Bresser, et al., 2020;Prasloski, Mädler, et al., 2012;Whittall & MacKay, 1989).One study reported that myelin quantification remains unaffected by the choice of N T 2 when it is sufficiently large, while a small N T 2 can result in substantial variations in MWF values ( Kumar et al., 2016( Kumar et al., , 2018)).In the current study, we chose N T 2 of 200 to maintain consistency with the number of single-T 2 values used Fig. 8. FLAIR images, T 2 maps, and MWF maps for three people with MS fitted using the data-driven approach.in creating the mcT 2 dictionary.The ability to generate a large mcT 2 dictionary in our study was made feasible by implementing a subsequent dilution process.This crucial step served to tailor the mcT 2 basis functions to match the physiological mcT 2 configurations found in the tissue and avoid erroneous combinations that could potentially arise due to noise.
The numerical results demonstrate that the datadriven algorithm can produce consistently accurate results, even in scenarios with substantial variations in the range of T 2 values, short T 2 fractions, and high B 1 + inhomogeneity levels.The numerical phantom was designed to mimic physiological MWF values, SNR levels, and matrix sizes matching in vivo scan settings.Accordingly, similar regularizations were applied to the numerical phantom and in vivo data.In this study, B 1 + estimation relied on MESE data rather than a separate B 1 + mapping scan.Recent study showed that acquiring an independent B 1 + map can improve the MWF results ( Mehdizadeh & Wilman, 2022).Nevertheless, the datadriven algorithm was able to produce relatively accurate B 1 + maps for all SNR values compared with the conventional approach, contributing to the accuracy of the final MWF maps.Further investigation is warranted to compare the reconstructed B 1 + maps using the data-driven approach with independently measured B 1 + maps.
Another important conclusion, demonstrated in the numeric simulation, was that the values of δ MV , and consequently ξ MV , should scale with SNR as higher noise levels require less stringent similarity criterion between simulated and experimental signals.Lastly, the stability maps presented in this study highlight the data-driven algorithm's robustness across a range of L 1 and Tikhonov regularization weights.Findings from experiments conducted on the physical phantom provide valuable insights into the performance of the data-driven algorithm in a controlled setting.The unique design of the phantom used in this study provided valuable ground truth reference.Despite certain model limitations, such as the absence of exchange and diffusion effects, this phantom offered a genuine benchmark to assess the data-driven algorithm's accuracy.The results demonstrate the data-driven algorithm's reliability in terms of both accuracy and precision, and across a wide range of MWF values.Scan-rescan experiments further emphasized the strengths of the data-driven approach, revealing high repeatability, exhibited as low interscan RMSD in comparison with conventional processing.It is worth noting that different regularization parameters were chosen for the physical and numerical phantom data.This adjustment was necessary due to the variations in scan parameters and physical properties of each phantom, such as voxel size and the distribution of T 2 values.These differences were carefully considered to ensure accurate evaluations.
The data-driven myelin mapping tool exhibited high reproducibility also in vivo.Repeatability tests demonstrated high correlation and small RMSD, which aligns with previous studies ( Levesque et al., 2010;Meyers et al., 2009).Examining MWF and T 2 maps versus the T 2weighted images from healthy subjects indicates that the intensity of qualitative T 2 -weighted images cannot be used as an accurate marker of myelin content, seeing as the intensity in qualitative images depends not only on myelin content, but also on other factors such as the total water content, macromolecular content, pathology, as well as on B 1 + inhomogeneities.Thus, from a radiological point of view, changes in T 2 weighting correspond more tightly with the tissue's single T 2 , rather than MWF values.Considering the comparison between the healthy subjects and people living with MS, data-driven MWF values were significantly different between healthy subjects and people with MS across all assayed NAWM ROIs, even after stringent correction for multiple comparisons.High AUC was furthermore achieved for all ROIs, indicating the potential of data-driven MWF values as a biomarker for MS.These findings align with previous studies ( Faizy et al., 2016;Laule & Moore, 2018) and indicate a consistent reduction in MWF values in people with MS in comparison with healthy subjects.In contrast to the phantom experiments, in vivo scans lack ground truth MWF values.The processing of in vivo data thus used regularization weights derived from the numerical phantom experiments, as they shared the same scan settings.
MWF values were estimated in this study also using a conventional approach, which differs from the data-driven approach by using a theoretical dictionary of single-T 2 sig-nals, rather than an mcT 2 dictionary of motifs which were learnt from the examined tissue.Corresponding MWF values were less accurate and precise for all assayed models.Numerical phantom results illustrated the conventional approach's inability to produce accurate values, exhibiting substantial mean absolute error across the entire range of regularization settings.Conventional analysis of the scanrescan data of the physical phantom (Fig. S6) showcased a wider spread of values, indicating increased sensitivity to random interscan signal variations and to noise vis-à-vis data-driven fitting.In the in vivo repeatability test (Fig. S8), the conventional approach exhibited similar correlation coefficient to the data-driven approach, but significantly higher MWF values across ROIs, similar to findings from the numerical simulations.Notably, this overestimation in healthy subjects persisted in people with MS as well (Fig. S9).Furthermore, in specific regions such as the frontal lobe and GCC, the conventional approach generated high AUC values albeit in an opposite direction to the natural progression of MS, showing higher MWF values for people with MS vs. controls (Fig. S10 and Table S3).
The primary cause of this overestimation can be attributed to the notably high Tikhonov regularization value, which was larger by two orders of magnitude compared with the data-driven approach.As described in Eq. ( 2), a higher Tikhonov regularization leads to a smoother spectrum, thereby necessitating the inclusion of more spectral components in the optimization process.Notably, conventional processing incorporated a broader range of T 2 values compared with the final set of selected motifs when using data-driven processing.It is important, however, to note that both data-driven and nondatadriven approaches produced identical T 2 maps, indicating that the center of mass of the spectra remained unchanged.The usage of a large Tikhonov regularization value causes selection of longer T 2 values, which was counteracted by increased energy in the short T 2 range.Utilizing a smaller regularization value (than the one which yielded the best results in the numerical simulations) would result in a significantly wider spread and lower quality in scan-rescan outcomes (result not shown).
Notwithstanding the successful results obtained using data-driven fitting, the current study has several limitations which should be considered.Firstly, the absence of ground truth in vivo limits one's ability to validate the calculated MWF values.While correlations with histology can provided some insights into microstructural features ( Laule et al., 2006), their applicability is constrained by postmortem tissue changes and the inherent variability in procedures such as fixation, slicing, and staining.Secondly, our proposed method, like other RNNLS-based approaches, operates under the assumption of a slow exchange regime, where the intercompartmental exchange is slow relative to the acquisition time (TE).Although previous reports suggested that exchange plays a minor role in MWF ( Kalantari et al., 2011), there are also claims that the exchange of water between microenvironments may bias MWF values ( Levesque & Pike, 2009), requiring to expand the tissue model to incorporate exchange, e.g., by using Bloch-McConnell equations ( Harkins et al., 2012).
Thirdly, in our experiments, maximal TE of 132 ms was used.Optimal estimation of the T 2 values of slow relaxing intra-/extracellular water pools would, however, require longer echo trains, which are not feasible in clinical settings due to the increase in specific absorption rate (SAR) caused by the addition of refocusing pulses.Furthermore, longer ETLs would limit the number of slices that can be acquired within one TR, once again requiring the use of longer TRs and longer scan time.Requiring a coverage of at least 9 cm within clinically feasible scan times thus forces the use of maximal TEs of 120-150 ms, thereby trading off some of the encoding quality of the longer T 2 components.Nevertheless, this constraint is less limiting when mapping MWF values seeing as we are mostly interested in the short T 2 component (T 2 < 40 ms), while trading off some of the accuracy in the estimation of the long T 2 components should have lower influence on the resulting MWF values.
It is important to note that MWF only serves as a proxy for myelin content, while reduction in measured MWF can, e.g., result from either lower myelin content (demyelination) or higher intra/extra water content (inflammation or edema).Previous studies have suggested that the total water content has a minimal effect on MWF, with reductions primarily attributed to demyelination ( Laule et al., 2004;Vavasour et al., 2021).Another factor that could influence the measured MWF is iron content, which is known to decrease T 2 relaxation time and may lead to overestimation of MWF ( Birkl et al., 2019).Contrasting findings, however, have been reported, as another study demonstrated lower MWF values in regions with higher iron content ( Khattar et al., 2021), suggesting that further investigation is needed to gain deeper insights into this cofactor.Additionally, incorporating other imaging contrasts alongside mcT 2 , such as diffusion ( Kolind et al., 2008;Rahmanzadeh et al., 2021) and magnetization transfer ( Schmierer et al., 2004), may provide additional valuable information and more accurate estimation of microstructural tissue changes.

CONCLUSIONS
This study introduces a new data-driven approach to mcT 2 analysis, which includes a new B 1 + correction procedure and incorporates entropy and pseudo-orthogonality regularizations in the simulated multicomponent signal model.
By drawing from concepts in statistics, the data-driven approach identifies global mcT 2 patterns in the WM, which are more likely to appear in the tissue, and are subsequently used to analyze the local signal in each voxel.This endows the resulting MWF values with high accuracy, precision, and robustness to noise when compared with the conventional RNNLS approach.The substantial difference in MWF values between healthy subjects and NAWM in people with MS indicates the potential of data-driven MWF values as a radiological biomarker for MS.
The ensuing values can be further extended to explore various aspects of the MS disease, including associated psychiatric conditions and cognitive impairment ( Curti et al., 2018), optic neuritis ( Reich et al., 2009), and remyelination ( Caverzasi et al., 2023).The data-driven approach itself can be also utilized for a broader range of applications involving multicompartment analysis, such as fat/ water separation ( Nassar et al., 2023), analysis of the prostate ( Storås et al., 2008), tumor characterization ( Nikiforaki et al., 2020), and generalized to other types of contrasts.

DATA AND CODE AVAILABILITY
The code used in this study will be shared publicly upon publication of the article.In vivo brain images presented in the study will be shared upon request, under approval of the local ethics committee and a data-transfer agreement.

Fig. 1 .
Fig. 1.Data-driven algorithm flowchart.Input for the algorithm is marked with dashed frames.(Step 1) Generating single-T 2 dictionary using the EMC algorithm.(Step 2) Creating mcT 2 dictionary by combining single-T 2 signals with different fractions.(Step 3) Signal correction is performed for compensating for transmit field inhomogeneities.(Steps 4, 5) Statistical correlation is computed between each dictionary element ( i and signal s j and added with entropy regularization to prevent overfitting.(Steps 6, 7) The result is then normalized and summed across all voxels to produce a global score.(Step 8) Select the basis elements with the highest scores while maximizing the orthogonality between motifs.

Fig. 2 .
Fig. 2. MWF mapping in a numerical phantom.Ground truth MWF and B 1 + maps are presented in the left panels.Right panels illustrate the performance of conventional and data-driven fitting approaches at different SNRs. Figure shows fitted MWF maps, absolute error maps, and reconstructed B 1 + profiles.

Fig. 3 .
Fig. 3. Stability of the conventional and data-driven MWF fitting approaches showing the mean absolute error (in %) as function of L 1 and Tikhonov regularization weights for numerical phantom across different SNRs.

Fig. 4 .
Fig. 4. Sagittal (A) and axial (B) illustrations of the unique multicompartment phantom design used in the study.(C) Correlation between data-driven and ground truth MWF values.Intratube variability (SD) of MWF values is shown as error bars for each tube.Black line denotes the best fit line.Correlation value of r = 0.99, RMSD value of 0.95%, and mean absolute error of 0.83 ± 0.51% were achieved between the estimated and ground truth MWF values.

Fig. 5 .
Fig. 5. Scan-rescan analysis and Bland-Altman plot for MWF values derived using the data-driven approach for two tubes containing (A, C) 7.1% and (B, D) 18.5%.

Fig. 6 .
Fig.6.T 2 -weighted images, T 2 maps, and MWF maps for three healthy subjects, generated using the data-driven fitting approach.

Fig. 7 .
Fig. 7. (A) Scan-rescan correlation analysis and (B) Bland-Altman plot for in vivo MWF values derived using the data-driven algorithm.

Table 1 .
List of optimal parameters for the conventional and data-driven approaches.