Improved myelin water imaging using B₁⁺ correction and data-driven global feature extraction: Application on people with MS

Myelin holds a pivotal role in the central nervous system (CNS), serving as a protective insulating sheath that envelops nerve fibers. This offers a crucial structural reinforcement and ensures the efficient transmission of electrical signals (Stadelmann et al., 2019). Consequently, the quantification of myelination levels is an invaluable biomarker with a wide range of applications. These include the study of heathy brain development, as well as investigation of neurodegenerative disorders such as Alzheimer’s disease (Dean et al., 2017), Parkinson’s disease (Dean et al., 2016), and multiple sclerosis (MS) (Laule et al., 2004; Laule & Moore, 2018; Loma & Heyman, 2011), where myelin mapping can enhance the understanding of these conditions, and facilitate their diagnosis, treatment, and disease progression.

Considering its significance, various MRI techniques have been developed for myelin water imaging (MWI), a reliable proxy for myelin content. The most prevalent among these methods are relaxometry-based techniques that rely on the difference between the rapid relaxation 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₁ and T₂ (mcDESPOT) (Deoni et al., 2013; Zhang et al., 2015), T₂ 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₂ relaxation times and avoid biases caused by field inhomogeneities, approaches based on spin echo protocols are favorable. Signal processing then involves multicomponent T₂ (mcT₂) analysis (Does, 2018; Whittall et al., 1997), where one maps the relative fraction of the fast-relaxing T₂ component (<40 ms) corresponding to water bound between myelin sheaths, and the longer T₂ 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₂ spectral content corresponding to myelin water to the total area under the T₂ spectrum.

The most efficient protocol for mcT₂ mapping in vivo is 2D (i.e., multislice) multiecho spin echo (MESE), offering high-sensitivity T₂ 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₂ 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₂ ‘motifs’) which best describe the examined tissue. These motifs are then used as basis functions for a regularized non-negative least square (RNNLS) mcT₂ 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₂ 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 (⁠$B1+$⁠) inhomogeneities, the use of entropy-based regularization constraint, and enforcing pseudo-orthogonality among mcT₂ 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.

For readers’ convenience, we provide a concise overview of the conventional, nondata-driven approach to mcT₂ 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₂ 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∈ℝETL$ is the experimentally acquired signal, $D∈ℝETL×NT2$ is a simulated dictionary of single T₂ EMC signals (⁠$di$⁠), $w∈ℝNT2$ is the relative fraction of each subvoxel water pool, ETL is the echo train length, and $NT2$ is the number of cellular components. Extraction of $w$ is typically done by solving an RNNLS minimization problem of the form:

where $λTikh,λL1≥0$ are Tikhonov and L₁ regularization terms. Tikhonov regularization is added to favor smoother solutions, while the L₁ regularization is added to promote sparse solutions. Due to the large space of possible mcT₂ 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.

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₂ 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₂ 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₂ EMC dictionary of signals, spanning $NT2$ logarithmically spaced T₂ values ranging from 10 to 800 ms. In step #2, a theoretical dictionary of all possible mcT₂ motifs is constructed by combining series of $NComp$ single-T₂ signals, each with a relative fraction $fn∈[0..1]$ using a fraction resolution of $Δf=1/Nfrac$⁠, where $Nfrac$ is the number of discrete values that $fn$ can assume. Each simulated mcT₂ dictionary element $d$ is expressed as

where $Ncomp$ 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₂ spectrum for each voxel, which will comprise a linear combination of mcT₂ motifs. An important substage of step #2 involves pruning nonphysiological configurations from the mcT₂ dictionary, including mcT₂ motifs that lack a short T₂ component (<40 ms); motifs whose short T₂ component fraction exceeds 30%; and motifs that correspond to a single-T₂ value that is outside the range of T₂ values that exist in the WM region. This single-T₂ value of an mcT₂ motif is calculated by fitting the motif’s signal decay curve to a dictionary of single-T₂ signals. Stage #3 consists of correcting for transmit field (B₁⁺) inhomogeneities, and is elaborated in Section 2.3.

In steps #4–5, statistical correlation is performed between each dictionary element $di$ and every voxel signal $sj$⁠, by calculating a “cost” $αi,j$⁠, defined as the L₂ 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):

such that

where $fi,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 (⁠$di$⁠) and voxel signal (⁠$sj$⁠). 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 $di$ and $sj$ 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 $di$ are summed across all voxels to yield a global score $Κi$ that expresses how well each motif matches the examined tissue:

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 pseudo-orthogonality 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.

After selecting the final set of motifs $D$⁠, the signal is modeled as a linear combination of mcT₂ signals:

where $s∈ℝETL$ is the experimental signal, $W∈ℝ|D|$ is the unknown vector of weights of the elements in $D$⁠, and $|D|$ is the number of elements in $D$⁠. The mcT₂ 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​∈ℝNT2×|D|$⁠, the final T₂ spectrum at each voxel will be

Finally, MWF values are calculated from the $T2$ spectrum of each voxel as the relative energy between 0 and 40 ms, and the energy of the entire spectrum.

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 $B1+$ profile and subsequently correct the experimental data. As each theoretical EMC decay curve $di$ depends on both $T2$ and $B1+$⁠, the same holds for each mcT₂ motif which can now be expressed as $d(B1+)$⁠. Thus, each motif $d$ transitions from a 1D to 2D by adding a $B1+$ dimension, discretized across $NB1+$ values in the range 80–120% (where 100% represents a fully homogeneous field).

The $B1+$ correction procedure involves three stages. First, the initial $B1+$ profile is calculated for each voxel j by finding the dictionary motif $d(B1+)$ that has the lowest L₂-norm difference to the experimental signal $sj$⁠. This is done using an exhaustive search over all dictionary elements and produces an initial solution for the transmit field map $B1,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, $B1,n+$ and $B1,n+1+$ denote the transmit field profile at iteration n and n + 1, respectively, and $B1+$ denotes the simulated range of $NB1+$ 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 $B1+$ value that best aligns with the data from voxel $sj$⁠; and second, a prior that imposes spatial smoothness within a 2D kernel using L₁ norm. The utilization of the L₁ 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, $Nk$ denotes all voxels within the 2D kernel surrounding voxel $sj$⁠, and $|Nk|$ is the number of voxels in $Nk$⁠. The iterative process is terminated either when there is no change in the $B1+$ value between iterations or when the number of iterations exceeds the predefined limit of $Niter=200$⁠. The resulting transmit field map is then denoted as $B1,opt+$⁠.

Finally, the $B1,opt+$ map is used to correct signal $sj$ in each voxel. Taking the mcT₂ motif with the optimal $B1+$ value $d(B1,opt+)$⁠, the corrected signal $sj,cor$ is calculated as

where $d$ denotes the “homogeneous” mcT₂ motif, corresponding to $B1+$ value of 100%.

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 (⁠$NComp=1,2,3$⁠), the myelin fractions, the relaxation times, and relative fractions of the intra-/extracellular water pools, and the $B1+$ field profile. Detailed description of simulated mcT₂ 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₂ fitting was performed according to the algorithm described in Figure 1, starting with generating an mcT₂ dictionary. The number of possible configurations (⁠$NmcT2$⁠) grows exponentially with $NT2$⁠, $Nfrac$⁠, and $NB1+$, 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 $NT2=200$⁠, $Δf=0.05$⁠, $NB1+=9$ (⁠$B1+=80:5:120%$⁠). The regularization weights $λTikh$⁠, $λL1$⁠, 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 $λL1$ values, for the five tissue types in the numeric phantom. Optimal reconstruction parameters for the conventional and data-driven approaches, i.e., those which yielded the highest accuracy, are detailed in the 3rd and 6th columns of Table 1, respectively.

An mcT₂ phantom was prepared using MnCl₂ solutions at concentrations of 0.11, 0.15, and 0.6 mM, producing T₂ 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₂ 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₂ (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.

MRI scans were conducted separately for each tube using a 2D MESE protocol with the following parameters: TE = 7.9, repetition time (TR) = 5000 ms, interecho spacing = 7.9 ms, ETL = 24, FOV = 750 x 750 mm², matrix size = 30 x 30, slice thickness = 3 mm, N_Slices = 10, bandwidth = 399 Hz/Px, and total scan time of T_acq = 2:45 min. A key aspect of these scans is the huge voxel size = 25 x 25 mm², ensuring that each tube was fully encapsulated within a single voxel. This unique setup enabled the generation of genuine, experimental mcT₂ signals, with known ground truth fractions of the short T₂ component. To generate a large number of mcT₂ 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₂ 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₂ 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.

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 (3933-17-SMC). MRI scans included a 2D MESE protocol with TE/TR = 12/5000 ms, interecho spacing = 12 ms, ETL = 11, FOV = 192 x 156 mm², matrix size = 192 x 156, slice thickness = 3 mm, N_slices = 32, acquisition bandwidth = 200 Hz/Px, GRAPPA acceleration factor = 2, and T_acq = 7:35 min. 3D T₁-weighted magnetization prepared two rapid gradient echo (MP2RAGE) data were collected for brain segmentations using TE/TR = 3.5/4000 ms, TI = 732 and 2220 ms, FOV = 224 x 168 mm², matrix size = 224 x 168, slice thickness = 1 mm, N_slices = 192, and T_acq = 6:24 min.

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² for $B1+$ correction. Single-T₂ maps were also calculated based on the mcT₂ spectrum within each voxel by (i) generating a single-T₂ signal using a linear combination of all T₂ components (see Eq. (11)), and (ii) performing single-T₂ 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.

A cohort of 29 individuals with relapsing–remitting MS (44.2 ± 11.8 y/o, 9 males) was scanned on a 3T scanner at Sheba Medical Center after providing informed consent and under Helsinki approval (6923-20-SMC). MRI scan included a 2D MESE protocol [TE/TR = 12/4600 ms, interecho spacing = 12 ms, ETL = 11, FOV = 192 x 220 mm², matrix size = 112 x 128, slice thickness = 3 mm, N_slices = 33, GRAPPA acceleration factor = 2, bandwidth = 200 Hz/Px, T_acq = 5:15 min]; 3D T₁w MPRAGE [TE/TR = 2.3/1800 ms, TI = 900 ms, FOV = 256 x 256 mm², matrix size = 256 x 256, slice thickness = 1 mm, N_slices = 176, T_acq = 4:35 min]; fluid attenuated inversion recovery (FLAIR) [TE/TR = 83/9000 ms, TI = 2500 ms, FOV = 240 x 195 mm², matrix size = 320 x 260, slice thickness = 2 mm, N_slices = 64, total scan time = 4:30 min]. MWF maps were generated using the same parameters as for the healthy subjects. T₂ maps were calculated from the mcT₂ spectrum at each voxel, following a procedure described above for the healthy subjects.

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.

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. $B1+$ 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 $B1+$ accurately and produced mean absolute errors of 6.0%, 5.1%, and 4.7%. At SNRs ≤ 100, higher deviations were observed for the $B1+$ maps using both approaches with mean absolute 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 $λL1$ regularization values and different SNRs. The data-driven approach exhibits consistently small mean absolute errors over a wide range of L₁ 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₂ 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₂ fitting is shown in Figure S6.

T₂-weighted images, T₂ maps, and MWF maps from three representative healthy subjects are presented in Figure 6. No correlation was observed between the MWF and T₂ values, i.e., higher MWF values do not necessarily indicate lower T₂ 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₂ fitting for healthy subjects is shown in Figures S7-S8.

Figure 8 shows FLAIR images, T₂ 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₂ 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₂ fitting is shown in Figure S10 and Table S3.

The significant ambiguity within mcT₂ 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₂ analysis of MESE data, with corrections for transmit field inhomogeneities. The method begins by identifying global mcT₂ motifs of the entire tissue, which are then used for localized signal analysis at each voxel, while accounting for local variations in the $B1+$ field. As an initial step, a comprehensive dictionary of mcT₂ 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₂ motifs significantly reduces the number of potential solutions, thereby alleviating the intrinsic ill-posed nature of mcT₂ 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 data-driven approach has been suggested (Piredda et al., 2022), where different regression models are trained using single relaxometry measurements (such as single T₁ and single T₂) 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₂ 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₂ dictionary does not impose any constraints on the final number of components in the T₂ spectrum, as there is no inherent limit on the number of mcT₂ motifs used to represent the signal in each voxel (see Eq. (11)).

The range of single-T₂ values used to filter nonphysiological mcT₂ 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₂ values in $D$ be sufficiently dense within the physiological range of T₂ values to construct an mcT₂ dictionary $D$ that will faithfully characterize the tissue. Various choices of $NT2$ 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 $NT2$ when it is sufficiently large, while a small $NT2$ can result in substantial variations in MWF values (Kumar et al., 2016, 2018). In the current study, we chose $NT2$ of 200 to maintain consistency with the number of single-T₂ values used in creating the mcT₂ dictionary. The ability to generate a large mcT₂ dictionary in our study was made feasible by implementing a subsequent dilution process. This crucial step served to tailor the mcT₂ basis functions to match the physiological mcT₂ configurations found in the tissue and avoid erroneous combinations that could potentially arise due to noise.

The numerical results demonstrate that the data-driven algorithm can produce consistently accurate results, even in scenarios with substantial variations in the range of T₂ values, short T₂ fractions, and high $B1+$ 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, $B1+$ estimation relied on MESE data rather than a separate $B1+$ mapping scan. Recent study showed that acquiring an independent $B1+$ map can improve the MWF results (Mehdizadeh & Wilman, 2022). Nevertheless, the data-driven algorithm was able to produce relatively accurate $B1+$ 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 $B1+$ maps using the data-driven approach with independently measured $B1+$ 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 $L1$ 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₂ 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₂ maps versus the T₂-weighted images from healthy subjects indicates that the intensity of qualitative T₂-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 $B1+$ inhomogeneities. Thus, from a radiological point of view, changes in T₂ weighting correspond more tightly with the tissue’s single T₂, 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₂ signals, rather than an mcT₂ 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 scan–rescan 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₂ 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 nondata-driven approaches produced identical T₂ 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₂ values, which was counteracted by increased energy in the short T₂ 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₂ 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₂ components. Nevertheless, this constraint is less limiting when mapping MWF values seeing as we are mostly interested in the short T₂ component (T₂ < 40 ms), while trading off some of the accuracy in the estimation of the long T₂ 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₂ 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₂, 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.

This study introduces a new data-driven approach to mcT₂ analysis, which includes a new B₁⁺ 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₂ 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.

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.

S.Z.: Methodology, Software, Validation, Formal analysis, Investigation, Writing—Original draft, Visualization; N.O.: Software, Investigation; D.R.: Investigation, Resources; N.S.: Investigation; T.B.-K.: Resources, Project administration; D.B.-A.R.: Resources, Validation; S.S.: Resources, Validation; C.H.: Resources, Validation; N.B.-E.: Conceptualization, Methodology, Software, Writing—Review & Editing, Visualization, Supervision.

This work was supported by Ichilov-TAU joint program for innovation in medical engineering (Alrov fund).

The authors declare no competing interests.

Supplementary material for this article is available with the online version here: https://doi.org/10.1162/imag_a_00254.

The minimization problem in Eq. (12) was solved using Matlab’s quadprog solver (Mathworks, version 2022b) for quadratic problems with linear constraints:

where H is semipositive definite matrix and q is a linear term. The optimization problem is cast into the quadratic programming solver by rearranging Eq. (12):

where $D$ is a dictionary of selected mcT₂ motifs, $W∈ℝ|D|$ is the unknown vector of weights of the elements in $D$⁠, $|D|$ is the number of elements in $D$⁠, and $λTikh,λL1≥0$ are Tikhonov and L₁ regularization weights. To keep only physical solutions, we constrain the elements in $W$ to be non-negative, i.e., $Wi≥0$ so $‖W1‖=∑i|Wi|=∑iWi$⁠. Rearranging the last expression yields

The term $sTs$ does not depend on $W$ and thus has no effect on the minimization and can be discarded:

where $I$ is the identity matrix and $1T$ is a vector of ones. Using the expression in Eq. (A.1), we can calculate that $H=DTD+2λTikI$ and $qT=λL11T​−STD$⁠.