Myelin water imaging from multi-echo T 2 MR relaxometry data using a joint sparsity constraint

Demyelination is the key pathological process in multiple sclerosis (MS). The extent of demyelination can be quantified with magnetic resonance imaging by assessing the myelin water fraction (MWF). However, long computation times and high noise sensitivity hinder the translation of MWF imaging to clinical practice. In this work, we introduce a more efficient and noise robust method to determine the MWF using a joint sparsity constraint and a pre-computed B1+-T2 dictionary. A single component analysis with this dictionary is used in an initial step to obtain a B1+ map. The T2 distribution is then determined from a reduced dictionary corresponding to the estimated B1+ map using a combination of a non-negativity and a joint sparsity constraint. The non-negativity constraint ensures that a feasible solution with non-negative contribution of each T2 component is obtained. The joint sparsity constraint restricts the T2 distribution to a small set of T2 relaxation times shared between all voxels and reduces the noise sensitivity. The applied Sparsity Promoting Iterative Joint NNLS (SPIJN) algorithm can be implemented efficiently, reducing the computation time by a factor of 50 compared to the commonly used regularized non-negative least squares algorithm. The proposed method was validated in simulations and in 8 healthy subjects with a 3D multi-echo gradient- and spin echo scan at 3T. In simulations, the absolute error in the MWF decreased from 0.031 to 0.013 compared to the regularized NNLS algorithm for SNR=250. The in vivo results were consistent with values reported in literature and improved MWF-quantification was obtained especially in the frontal white matter. The maximum standard deviation in mean MWF in different regions of interest between subjects was smaller for the proposed method (0.0193) compared to the regularized NNLS algorithm (0.0266). In conclusion, the proposed method for MWF estimation is less computationally expensive and less susceptible to noise compared to state of the art methods. These improvements might be an important step towards clinical translation of MWF measurements.


Introduction
Myelination is a crucial aspect of brain development and is essential for the functioning of the nervous system. Demyelination, on the other hand, is a pathological process that plays an important role in certain diseases such as multiple sclerosis (Laule et al., 2004;Faizy et al., 2016). Accurate measurement of myelin content has the potential to increase our insights into several disease processes. Magnetic resonance imaging (MRI) enables imaging of features related to (de-)myelination in vivo. Methods to do so include ultra-short echo time imaging, diffusion tensor imaging, magnetization transfer imaging and multi-echo T 2 (MET 2 ) or T * 2 relaxometry methods (MacKay and Laule, 2016;Does, 2017). In MET 2 relaxometry (Whittall and MacKay, 1989;Poon and Henkelman, 1992;Mackay et al., 1994), a T 2 distribution is determined from a multi-echo spin-echo (MESE) acquisition. The analysis of this distribution is normally limited to the white matter, in which the short T 2 relaxation times (10-40 ms) are considered as myelin water (MW), intermediate T 2 relaxation times (40-200 ms) as intra-and extracellular water (IECW) and longer T 2 relaxation times (>1s) as free water (MacKay and Laule, 2016). The myelin water fraction (MWF) is calculated as the ratio between the signal contribution from MW to the total sum of signal contributions. It was shown that the method results in reproducible MWF maps, but these maps are dependent on methodological variability (Levesque et al., 2010;Meyers et al., 2013).
A Gradient-and Spin-Echo (GRASE) acquisition pattern (Oshio and Feinberg, 1991) was introduced in the field of ME T 2 relaxometry to obtain whole brain images with shorter acquisition times compared to a regular MESE acquisition (Prasloski et al., 2012b). The outcomes of regular MESE and multi-echo GRASE imaging are known to be highly similar (Kumar et al., 2018;Drenthen et al., 2019a).
However, MESE and multi-echo GRASE signals are sensitive to inhomogeneities in the B 1 transmit field. In particular, these inhomogeneities can cause refocusing pulses with flip angles very different from 180 ∘ , leading to secondary and stimulated echos (Hennig, 1991). As a result, deviations from mere exponential decays can occur. To account for B 1 inhomogeneity effects, the signal for each T 2 component may be calculated based on the extended phase graph (EPG) formalism (Hennig, 1988) using the corrected flip angle. Such flip angle inhomogeneity (FAI) correction is especially important at higher field strengths (B 0 ! 3T), since in that case the B 1 field becomes less homogeneous.
In addition, the inverse problem of computing a T 2 distribution from MET 2 data is highly underdetermined and therefore very sensitive to noise (Graham et al., 1996). Therefore, it is crucial to acquire data with a high signal to noise ratio (SNR) and apply regularization in order to obtain a stable solution. A common approach (Whittall and MacKay, 1989) is to solve the problem for every voxel independently assuming that the T 2 distribution is smooth. This smoothness is enforced by including the first-order derivative of the T 2 distribution as a penalty term in the objective function. The resulting problem can be solved by the regularized non-negative least squares (regNNLS) algorithm (Lawson and Hanson, 1974). To include FAI correction, Prasloski et al. (2012a) proposed an EPG approach using regNNLS fitting, which was solved for different FAI. This makes it possible to find the optimal FAI value to calculate the final T 2 distribution in a voxel-wise manner. However, including this correction increases the computational complexity leading to very long processing times. Yoo et al. (2015) were able to reduce the computation time by a factor 4 using CPU and GPU parallelization leading to 10 min computation time for a dataset with matrix size 256Â 256 Â 7 and 32 echos. Although it has been shown that this method gives reproducible results (Lee et al., 2018;Drenthen et al., 2019b), it remains sensitive to noise, leading to relatively large coefficients of variation. As such, the method can benefit from improved regularization or other techniques that better control the noise amplification.
Several other regularization approaches were proposed, either leading to simplifications of the signal model or to increased computation times. In Hwang and Du (2009), it was proposed to include additional (2D) spatial regularization in the regNNLS algorithm, leading to smoother MWF maps. Similarly, 2D spatially regularized MWF mapping with B 1 inhomogeneity correction was proposed by Kumar et al. (2016) and later extended to include spatial smoothness of the FAI map and 3D spatial regularization (Kumar et al., 2018). Due to the complexity of the problem, the computation time was reported to take approximately 15-16 h for MET 2 data of size 80 Â 80 Â 64 with 32 echoes. While the above mentioned methods allow for a large number of T 2 components (up to 50 or 100), Stanisz and Henkelman (1998); Raj et al. (2014); Hajj et al. (2019) proposed a two or three compartment model, assuming that the distribution can be described by Gaussian peaks representing MW, intra-and extracellular water (IEW) and free water, respectively. Constraining these components to a predefined T 2 range combined with the restriction of having a fixed, small number of components drastically reduces the flexibility of the model, but yields improved noise robustness. In a similar way, Akhondi-Asl et al. (2014) used an inverse Gaussian distribution for a three compartment model, which prevents a distribution with a tail reaching to negative T 2 values.
Recently, several methods were introduced that used a T 2 distribution consisting of delta peaks. Bj€ ork et al. (2016) showed that the distribution does not necessarily need to be smooth, since a Gaussian distribution or combination of delta peaks essentially lead to the same measurement. For exponential signals, this makes it possible to estimate the MWF in a parameter-free manner, without a pre-defined T 2 grid through a system identification approach (EASI-SM algorithm) (Stoica and Babu, 2013). However, since this algorithm is specifically designed for exponential signals, it does not allow the correction for FAI using EPG simulations. Drenthen et al. (2019b) recently proposed to use orthogonal matching pursuit (OMP) instead of the regNNLS algorithm. They applied the non-negative OMP algorithm proposed by Yaghoobi et al. (2015) and Nguyen et al. (2017), which implicitly includes a non-negativity constraint. The NNLS and NNOMP algorithm show great similarities, but Drenthen et al. demonstrated that applying temporal regularization in the NNLS algorithm leads to a bias in the estimated MWF. Recently, Does et al. (2019) proposed a method based on principle component analysis to distinguish the components contributing to the signal from those characterizing the noise. In this manner the method provides a way to pre-process noisy relaxometry data.
In most of these studies, the proposed algorithms were compared to the state of the art regNNLS algorithm. In these comparisons, several methods showed a higher MWF in the sub-cortical white matter and major white matter tract regions, which was confirmed by the signal in T 2 weighted scans.
Most recently, dictionary-based methods have gained increased interest for quantitative MR parameter mapping. Popular examples include applying a dictionary as a signal representation for compressed sensing image reconstruction (Doneva et al., 2010;Li et al., 2012), a grid search for fast parameter estimation (Marques et al., 2010;Barral et al., 2010;Ben-Eliezer et al., 2015), and MR Fingerprinting (MRF) in which multiple tissue parameters and system parameters are estimated simultaneously (Ma et al., 2013). In these methods, a pre-computed dictionary containing simulated signal evolutions is used for tissue and/or system parameter mapping. Specifically, the inner product is applied to identify the best matching dictionary atom and in effect the corresponding parameters from a measured signal. The dictionary depends on the pulse sequence and the expected range of tissue and system parameters that need to be estimated. This dictionary is computed once for a given pulse sequence and can be reused for all subsequent acquisitions.
Recently, several algorithms were proposed to perform a multicomponent analysis of MR Fingerprinting data (McGivney et al., 2017;Tang et al., 2018), where components are distinguished based on T 1 and T 2 values. Similarly to multi-component T 2 approaches, these methods perform the multi-component MRF (MC-MRF) analysis for each voxel separately applying a sparsity constraint to limit the number of components per voxel. Following on these methods, we recently proposed a new method for MC-MRF based on the NNLS algorithm that applies a spatial joint-sparsity constraint leading to a small number of components across the region of interest (Nagtegaal et al., 2020). This additional constraint enables further noise resilience of the estimated component weights and implementation in a computationally efficient algorithm. Consequently, it leads to significantly reduced computation time compared to the Bayesian and reweighted ℓ 1 -norm approaches.
In this work, we propose a new multi-component approach to MWF mapping, based on our previously proposed algorithm for MC-MRF, which is extended to include FAI correction. The FAI map is initiated by performing a voxel-by-voxel, dictionary based, single-component parameter estimation. Subsequently, a multi-component analysis is performed with an algorithm combining a joint-sparsity constraint with nonnegativity, in which the correction for FAI effects is included. We assume that the T 2 distribution is sparse and all voxels within a region of interest share the same T 2 components. This is a crucial difference with the common assumption that the T 2 distribution is temporally, and (optionally) spatially smooth. We hypothesize that our approach will reduce noise amplification in the MWF maps and allow for a computationally efficient algorithm. The proposed method is evaluated in numerical simulations and in vivo measurements and compared to the regNNLS algorithm.

Data model
The multi-component signal x j 2 R M measured in voxel j at M time points for a MESE sequence is modeled as where Sða j ; T 2 Þ 2 R M is the signal for relaxation time T 2 at FAI value a j , c j ð ÁÞ is the T 2 distribution, and e j a Gaussian noise vector. The value of c j ðT 2 Þ can be considered the signal contribution of a tissue with a certain T 2 time in voxel j.
The MESE signal decay Sða j ; T 2 Þ for non-ideal refocusing pulse flip angles can be calculated using the EPG formalism (Hennig, 1988). The applied sequence consists of an a j Á 90 ∘ pulse followed by M a j Á 180 ∘ pulses.
The integral of Equation (1) might be discretized by taking N T2 T 2 -values and N FAI FAI values assuming that there are N ¼ N T2 Á N FAI possible signal realizations. These signals will be stored in a matrix D 2 R MÂN , to which we will refer as the dictionary. A subdictionary containing the signals for a specific FAI value a is indicated as D a 2 R MÂNT 2 . In this matrix the rows correspond to the simulated signal and each column to a particular T 2 time.
Assuming that there is only a single FAI-value in voxel j, Equation (1) can be written as a linear combination of N T2 signals with weights c j , corresponding to the discretized T 2 distribution values and the FAI value a j for the given voxel: The weights of the T 2 distribution are assumed to be non-negative. Given a measured signal x j 2 R M the weights of the T 2 components can then be estimated by solving the non-negative least squares problem: Conventional methods for solving this minimization problem assume that the vector c is either smooth or sparse. Very recently, we introduced a different approach by imposing a joint sparsity constraint (Nagtegaal et al., 2020). Our premise was that the measured signals in a region of interest (ROI) could be described by a small set of T 2 relaxation times, common for all voxels in the ROI. To formalize this, let C ¼ ½c 1 …c J be the N T2 Â J sized matrix containing the contributions for all J voxels. Furthermore, c i is taken to represent a row from this matrix, corresponding to the contributions of a particular T 2 signal to all voxels, so that at the same time C ¼ ½c 1T …c NT 2 T T .
This leads to the joint sparsity minimization problem: The term P NT 2 i¼1 kc i k 0 counts the number of used components, but does not restrict this to a prescribed maximum. The joint sparsity problem including the non-negativity constraint can be solved with different algorithms. We applied the Sparsity Promoting Iterative Joint NNLS (SPIJN) algorithm (Nagtegaal et al., 2020) for MC-MRF involving highly-coherent signals. This approach enabled a higher noise robustness compared to voxel-by-voxel methods and lead to easy interpretable results because of the small number of components. In this work we extend the SPIJN algorithm to include FAI correction and investigate its application to multi-component T 2 analysis.
Algorithm 1. The Sparsity Promoting Iterative Joint NNLS (SPIJN) algorithm to perform a multi-component analysis for MET 2 data with correction for flip angle inhomogeneities.

Fitting procedure
We propose a two-step approach to perform multi-component T 2 analysis of MET 2 data: (1) a FAI map is computed assuming that the measured signal is dominated by a main component, and thus can be modeled as a single component in each voxel; (2) a multi-component T 2 analysis is performed using the estimated FAI map and applying the joint sparsity constraint as stated in Eq. (4). The algorithm is schematically described in Algorithm 1.
Using the EPG formalism a dictionary D was computed containing combinations of FAI and T 2 -values. A fixed T 1 ¼ 1 s was used in the EPG simulations similarly to (Prasloski et al., 2012a;Kumar et al., 2018). FAIs were simulated as a multiplicative factor modifying the prescribed flip angle. The FAI values are modeled as a ¼ αeffective αintended , for a CPMG sequence the effective signals are symmetric around a ¼ 1. The modeled FAIs ranged from 0.75 to 1 in 140 linear steps, while T 2 relaxation times were chosen on a logarithmic scale from 10 ms to 5 s with 141 steps. The total computation time for the dictionary was 81 s.
The dictionary was first used to perform a single component matching. For each voxel, the inner product between measured signal and (normalized) dictionary signals was used to determine the FAI-T 2 combination that best described the measured signal.
Subsequently, the same dictionary was applied in the SPIJN algorithm for multi-component estimation and while doing so the FAI was restricted by the value obtained through single component matching. The SPIJN algorithm was based on the NNLS algorithm and used an iterative reweighting scheme to couple the non-negative sparse solutions of the different voxels. By applying this reweighting the solution converged to a jointly sparse solution.
The proposed method was implemented in Python. More details on the used convergence thresholds, the used reweighting and regularization can be found in (Nagtegaal et al., 2020). The regularized NNLS algorithm including FAI correction (Prasloski et al., 2012a) was used as a reference method. 101 T 2 values logarithmically spaced from 10 ms to 5 s and a fixed T 1 of 1 s were used. The regNNLS computations were performed with MATLAB 2018b (The MathWorks Inc; Natick, Massachusetts, USA).

Numerical simulation experiments
Two numerical experiments were performed to analyze the proposed method. Both experiments were performed for M ¼ 48 echoes with a spacing of ΔTE ¼ 10 ms (first echo at 10 ms). These settings were also applied in the in vivo experiment (see below).
First, the behavior of the proposed FAI estimation was analyzed and compared to FAI estimation with the regNNLS algorithm. Simulations were performed based on a signal composed as the weighted sum of two components: a short water (SW) relaxation component, with T 2 ¼ 20 ms, and a long water (LW) relaxation component with T 2 times from 25 ms to 3 s with 41 steps on a logarithmic scale. The SW fraction (SWF) ranged from 0 to 1 with step size 0.05 (while LWF ¼ 1 -SWF). The FAI level varied between 0.75 and 1 in 5 steps. Furthermore, 100 real valued Gaussian noise realizations were added to each simulated signal, to yield 100 noisy signal versions at each setting. The Gaussian noise had a standard deviation defined as s 0 =SNR, in which s 0 was the signal intensity of the first echo. A fixed SNR of 250 was used, which was comparable to the SNR of the in vivo experiment (see below). The absolute value of the noisy signal was then analyzed with the proposed FAI estimation. The mean residual signal error as well as the mean FAI error were calculated for each parameter combination (SWF/FAI/LW-T 2 ).
The mean residual signal error was computed as 1 100 P 100 i¼1 ks À d i k 2 = ksk 2 where s is the ground truth noise free signal and d i is the matched dictionary signal for noise realization i. Second, the precision and accuracy of the calculated MWF with the proposed method was compared to outcomes of NNLS and regNNLS algorithms. Therefore, an image was simulated consisting of 100 Â 100 pixels with a mixture of three components. These components had T 2 relaxation times of 20 ms, 70 ms and 1000 ms roughly corresponding to MW, IECW and FW, respectively. The map with the signal fractions (summing to one) of the components is shown in Fig. 1.
The simulations were performed with two different, constant FAI levels of 1.0 (i.e. no offset in B 1 field) and 0.9. The simulations were performed with different noise realizations, applying SNRs of 500, 250 and 100.
The MWF was computed with the unregularized NNLS and SPIJN algorithms using the FAI map estimated in an initial single component analysis step as described above. For the SPIJN algorithm λ ¼ 0:02 was used. These MWF maps were compared to the MWF from the regNNLS algorithm and to the ground truth. Relaxation times below 40 ms were considered to correspond to the MW component.

In vivo imaging experiments
To study the practical feasibility of the proposed method, brain scans were performed in 8 healthy subjects. Informed consent was obtained from all subjects. Imaging was performed on a 3.0 T Ingenia scanner (Philips, Best, The Netherlands) based on a 3D multi echo GRASE acquisition scheme using a 13 channel head coil. The sequence parameters were: 48 echoes with echo spacing of ΔTE ¼ 10 ms; EPI factor of 3; field of view 240 Â 205 Â 72 mm 3 ; voxel sizes 1:25 Â 1:25 Â 8 mm 3 ; repetition time TR ¼ 1:2s, resulting in a total acquisition time of 6 min and 14 s.
A multi-component analysis of the acquired brain data was performed using the proposed method. The proposed method was applied to all slices simultaneously. Skull and air were masked based on their signal intensity and other tissues not connected to the brain were also discarded as such. Additionally, the regNNLS algorithm was applied for myelin water mapping. For the SPIJN algorithm, the regularization parameter λ was set to 30 for the in vivo data (see Supplementary material Figure S1 for a range of λ values). The part of the T 2 distribution with a relaxation time shorter than a preset threshold valueT 2 was attributed to MW. The MWF was calculated as MWF ¼ P Experimentally, we studied two different thresholds:T 2 ¼ 30 ms and T 2 ¼ 40 ms. Regions of interest (ROIs) were manually annotated: splenium, thalamus and genu of the corpus callosum and parts of the frontal, occipital and temporal white matter lobes. The SNR of the in vivo data was calculated from the results of the regNNLS algorithm as the ratio between the signal intensity in the white matter structures in the first echo and the standard deviation of the residual.
For each method, subject and ROI the mean MWF and coefficient of variation (CoV) were calculated. A paired two-sample t-test was used to determine the significance of the differences in MWF values and CoV per region.  single component signal approximation error for varying FAI, SWF and LW-T 2 relaxation times (averaged over 100 noise realizations) for the SPIJN and regNNLS algorithm. The signal approximation error is not shown for the regNNLS, since it was negligible (max. 5 Á 10 À4 ). The middle row demonstrates that the error in the proposed method is largest for a combination of a SWF between approximately 0.5 and 0.9, large LW-T 2 relaxation times and more severe flip angle inhomogeneity level (FAI toward 0.75). Furthermore, the bottom images show that when the SWF and LW-T 2 increase, the signal estimation error increases as well because the signal model becomes less accurate. A LW-T 2 < 160 ms and SWF 0:2 is a realistic range for a mixture of MW and IECW, see e.g. (MacKay and Laule, 2016). In this range, demarcated by the dotted line in the figure, the maximum mean absolute error in the estimated FAI is 0.0266.

Numerical simulation experiments
Compared to the regNNLS algorithm (first row) the proposed method shows similar errors in the range of interest, but not for more extreme combinations. The overall mean FAI error using the regNNLS algorithm was lower (0.0158) compared to the proposed method (0.0460). However, the maximum error in the realistic range was 0.0301 with the regNNLS method, which is higher than the error of the proposed method (0.0266). Also, the mean error in this range was higher: 0.0146 compared to 0.0133, for regNNLS and proposed approach, respectively. In Fig. 3 the absolute mean MWF error maps (over 100 noise realizations) are shown for a FAI values of 1 and 0.9, summarizing the second MWF estimation experiment. Above each map, the root mean square error (RMSE) is indicated. Estimates were calculated with the NNLS, SPIJN, and regNNLS algorithms, at SNR values of 500, 250 and 100.
The SPIJN algorithm resulted in lower MWF RMSE for all SNR values compared to the NNLS and regNNLS algorithms. Observe that especially for FAI ¼ 1 the RMSE were markedly lower. Particularly for SNR ¼ 250 the SPIJN algorithm achieved a 58% lower error compared to the regNNLS algorithm.

In vivo imaging experiments
The mean SNR of the in vivo MET 2 data sets, as defined in Section 2.4, was 320 (maximum: 436; minimum: 254). Application of the SPIJN algorithm led to 5 to 7 components. The distribution of T 2 -values across identified components in the subjects is shown in Fig. 4. In general, it can be seen that components are matched to the lower and upper bound of T 2 values in the dictionary (10 ms and 5 s respectively). In between these bounds, components are matched to T 2 values around 35, 75 and 250 ms. The fraction maps for the different components for a representative subject (subject 2) are shown in Fig. 5.   Fig. 3. Absolute error maps of estimated MWF with the NNLS, SPIJN and regNNLS (Prasloski et al., 2012a) algorithms for FAI values of 0.9 and 1 in simulations with different SNR. Underlying ground truth fraction maps are shown in Fig. 1. The maps of the NNLS and SPIJN algorithms were computed using FAI estimations from the proposed single component dictionary matching. The root mean square error (RMSE) for each map is indicated on top. Fig. 4. T 2 value distribution of components identified by the SPIJN algorithm for 8 subjects. The typical threshold value for MW detection (T 2 ¼ 40 ms) is marked by a red line. The distance between the light grey grid lines reflects the dictionary step size. The size of the dots shows the relative abundance of the different components. Fig. 6 shows representative T 2 weighted images, estimated FAI maps and MWF maps (T 2 ¼ 40 ms) for the same subject. FAI and MWF maps were obtained with the SPIJN algorithm and the reference regNNLS algorithm. The FAI maps of the two methods show good agreement. The main difference between the regNNLS and SPIJN MWF maps can be seen in the frontal white matter (red circle). Here, the T 2 weighted images show evidence that there is myelin present, signified by the white-grey matter contrast. However, the regNNLS algorithm estimates very low MWF in these regions. Fig. 5 shows that the T 2 component of 36.2 ms is in particular responsible for the higher MWF detected by the SPIJN algorithm.
The mean MWF values across the ROIs and the subjects for the regNNLS and SPIJN withT 2 ¼ 30 ms and 40 ms are summarized in Table 1. Furthermore, CoVs are collated in Table 2. For reference, literature values as measured with different algorithms are shown in Table 3. Observe that forT 2 ¼ 30 ms differences between the methods are not significant. However, they are significant in most structures forT 2 ¼ 40 ms. Only for the splenium the outcomes of the two methods are not significantly different at this threshold. Simultaneously, note that for T 2 ¼ 40 ms the CoVs are significantly smaller for the proposed method compared to the regNNLS method.
The distribution of the differences between the two methods across the regions forT 2 ¼ 30 ms are plotted in Fig. 7. In most subjects the differences are not significant; only in subjects 7 and 8 significant differences were found at this setting. Notably, these two subjects yield components with T 2 ¼ 25 ms (see also Fig. 4), which adds to the MWF at T 2 ¼ 30 ms. The average computation time per slice was 1.19 s for the single component matching and 7.00 s for the multi component matching with the SPIJN algorithm. The average computation time per slice for the regNNLS was 408 s.

Discussion
In this study the SPIJN algorithm was introduced as a new, fast method to determine the MWF from MET 2 relaxometry data through a multi-component analysis with a correction for flip angle inhomogeneity. The method was compared to the NNLS and state of the art regNNLS algorithms in numerical simulations and on data from 8 subjects acquired at 3 T.
The first simulation experiment showed that a FAI map can be accurately estimated, especially in a realistic range of IECW relaxation times and MWFs. The FAI-mapping method is proposed as a computationally efficient technique and relies on the presence of a dominant tissue or small differences in T 2 between different components. The regNNLS approach (Prasloski et al., 2012a) provides a more generally applicable method, but is computationally much more expensive. The mean error of the estimated FAI increases (order of 20%) with increasing FAI, MWF (> 0:4) and increased T 2 of the non-MW component (> 400 ms). This could result in a biased FAI estimation e.g. in partial volume voxels with mainly myelin and cerebrospinal fluid. However, this setting represents a non-realistic configuration in myelinated tissue.
The second simulation experiment as shown in Fig. 3 demonstrated that the SPIJN algorithm yields lower RMSE for MWF estimation than the NNLS and regNNLS algorithms. The shown simulations used a signal model with a sparse distribution of T 2 values, possibly favoring SPIJN over the regNNLS algorithm, which assumes a smooth distribution. As such, the difference between SPIJN and NNLS, in which a sparse distribution is the only assumption, particularly shows the benefit of the joint sparsity constraint. The assumption of a jointly sparse T 2 distribution makes it possible to benefit from shared information in different regions and large numbers of voxels. As a consequence, improved noise resilience could be observed, especially at SNR ¼ 100. In-vivo measurements this may be important for balancing the trade-off between resolution and SNR. Reduced SNR caused by higher resolutions are expected to have smaller effects on accuracy with the proposed SPIJN algorithm compared to voxel-wise methods.
The number of resulting components enforced by the joint sparsity constraint may be difficult to predict a priori due to the complexity of tissue microstructure and potential natural variation. The number of retained components is influenced by the T 2 values of tissues present in the region of interest, the sensitivity of the used MET 2 sequence to these tissue parameters and the used level of regularization. As shown in Fig. 4 the in vivo experiments yield small variations in the number of estimated components and associated T 2 relaxation times, which is probably due to natural diversity.
In general the T 2 values of the main groups of components are in agreement with the myelin water, intra-extracellular water and free water, as observed on a voxel basis by e.g. Whittall et al. (1997). Three components with T 2 < 40 ms were identified in 4/8 subjects, while only two such components were obtained in the other 4/8 subjects. Fig. 5 illustrated that the additional component (T 2 ¼ 20:3 ms in this case) typically was associated with small signal fractions. In addition to these MW components, all subjects yielded a component with T 2 around 70 ms, which forms the main component in most voxels and is attributed to intra and/or extra cellular water. The identification of two further components with longer relaxation times is consistent with our earlier findings with MC-MRF (Nagtegaal et al., 2020). The longest T 2 component (T 2 ¼ 5s) can be contributed to free water and is mainly present in the locations where free water is expected. The component around T 2 ¼ 250 ms is rather small (typically 10%), but is in confirmation with earlier observations (Whittall et al., 1997;Laule et al., 2007b) and was interpreted as extra-axonal water (Does, 2017). A dominant component with comparable relaxation times was observed in patients with MS (Laule et al., 2007a). This component is not often reported in healthy subjects. The appearance of this component could be caused by increased sensitivity of the proposed algorithm to components with a low presence and thus contain useful information. It could also be caused by the proposed analysis method and then it has to be considered as an artifact. The information contained by this component would be of interest in further clinical studies or in an analysis of already acquired data.
Underlying biological principles could also cause a slight variation of T 2 values of the same component. This would conflict with the assumption of group sparsity, which determines a small basis to represent all measured signals. Since small differences in T 2 only cause a small difference in the exponential decay the effect of a slight variation on the T 2 -components are minor. This was further assessed in an experiment as described in Supplementary material Figure S2.
The focus of the rest of the paper is on the accuracy with which a MWF map could be obtained from components reflecting low T 2 times. The values obtained with the SPIJN algorithm atT 2 ¼ 40 ms are in agreement with the literature values as given in Table 3. Only the mean MWF in the genu of the corpus callosum shows a large difference with the values as reported by Kumar et al. (2018). On the other hand, they are in line with the results as reported by Drenthen et al. (2019b). The MWF Fig. 6. T 2 -weighted images, FAI and MWF maps for the regNNLS and SPIJN algorithm for a representative subject. T 2 40 ms was considered to correspond to MW. Frontal white matter regions with increased white-grey matter contrast are annotated with red ellipses. In slice 4 an example of the manually annotated ROIs is shown.
values for the regNNLS algorithm are consistent with the values as reported in (Prasloski et al., 2012b). The MWF maps obtained by the SPIJN algorithm showed higher fractions compared to those computed with the regNNLS approach. This difference might be due to the smoothness constraint on the T 2 distribution imposed by regNNLS. In effect, myelin water with a T 2 around 35 ms could be smoothed into the larger IECW pool, so that it is not identified as MW.
Notably, the shown T 2 weighted images indicated that higher levels of myelin might be expected than reflected in the regNNLS maps. Moreover, the MWF maps obtained by the SPIJN algorithm are very similar to MWF maps presented by Kumar et al. (2018). The main difference compared to the regNNLS algorithm can be seen in the frontal white matter. A possible lack of sensitivity in this region was reported by Wiggermann et al. (2019) as well, where it was hypothesized that this insensitivity was caused by the increased flip angle inhomogeneities in this region. The estimation of FAI presented by us is slightly different and combined with the improved noise resilience this could result in an improved MWF estimation in these regions.
The work by Kumar et al. (2018) also reported MWF values in several regions of interest that are comparable to those found with the SPIJN algorithm. Unfortunately, only a small number of studies reported MWFs for different structures, making a reliable comparison difficult. Possibly the limited number of studies reporting this is because of large natural variation in subjects. When values were reported (see e.g. (Prasloski et al., 2012b;Kumar et al., 2018;Drenthen et al., 2019b)), this was done as part of the introduction of new acquisition methods or algorithms, using only a small number of subjects in a similar way as was done here.
The regularization parameter in the SPIJN algorithm (λ, see Algorithm 1) was experimentally determined, but showed to be robust across different scans from the same scanner, once the same λ was used for processing of all scans after the value was set. Variations in the regularization parameter (from 20 to 60; observe that 30 was used in the results) only had very small effects on the resulting MWF maps (less than 1%, see Supplementary material Figure S1). In the simulations a different regularization parameter was used compared to the in vivo data. This was necessary because the simulations were based on a different number of tissue types and associated T 2 times.
The regularization level mainly determines the number of components and only indirectly the fractions of the different components. This makes the SPIJN method less sensitive to the exact regularization value. Essentially, setting lambda smaller will lead to a reduced fit error, but simultaneously to a larger number of components. Consequently, the regularization parameter might for example be automatically selected by requiring that a certain minimal fit error is achieved. However, performing such optimization for every scan will go at the expense of increased computation time.
We propose to use a jointly reweighted NNLS scheme to approximate the multi-component problem from (4). Other optimization schemes could be considered e.g. group LASSO (Yuan and Lin, 2006) as well. However, we experienced that the highly coherent T 2 dictionary atoms strongly affect the convergence properties of some of these methods. This confirms findings by others that FISTA or LASSO based methods are not always suitable choices (Wipf and Nagarajan, 2010;Tang et al., 2018).
Very recently two papers were published on the use of deep learning networks for the calculation of MW fractions based on the regNNLS algorithm (Lee et al., 2020;Liu et al., 2020). These networks are trained based on the regNNLS algorithm and are therefore able to reproduce these maps very well and therefore show the same level of noise sensitivity as the regNNLS algorithm. Both papers applied this in a Significant differences between the two methods (p < 0:05) are indicated by *. Significant differences (p < 0:05) between the two methods are indicated by *. Table 3 MWF values across several brain structures as reported in literature. All methods usedT 2 ¼ 40 ms. Studies with * were performed at 1.5 T, other studies were performed at 3 T as the results shows in this paper.
Brain structure (Whittall et al., 1997 Fig. 7. Box plots of difference between the MWF obtained for the SPIJN and regNNLS algorithms in 7 anatomical structures across the eight subjects. T 2 30 ms was considered to correspond to MW.
voxel-by-voxel manner, where the here proposed method aims to improve the model through the addition of the joint sparsity assumption. The proposed method was here applied on ME T 2 relaxometry data. The use of a joint-sparsity constraint could also be beneficial for other T 1 , T 2 or T * 2 relaxometry methods, as long as a similar multi-component model is applicable. The here proposed method was only evaluated in simulations and healthy subjects with GRASE acquisitions, focusing on the differences in the normally expected signals and MWF maps. We consider the evaluation of the method in patients a very important topic for future research.

Conclusion
The SPIJN algorithm facilitates estimation of MW fractions through a joint sparsity promoting fit of multiple T 2 components to T 2 relaxometry data. The method yielded enhanced accuracy in simulations compared to the state-of-art regularized NNLS algorithm. Furthermore, the MWF maps from healthy subjects showed visual improvements over the regNNLS approach. The method was also 50 times faster than the regNNLS algorithm: the average computation time per slice was 8.19 s on a standard desktop PC.
The faster computation of MWF maps combined with improved accuracy can help to increase our insights into (de)myelination and enables reconstruction of MWF maps directly after data acquisition.