Synopsis

Keywords: Quantitative Imaging, Multi-Contrast, Myelin

Myelin water fraction (MWF) mapping using MRI has enabled researchers to directly examine myelination and demyelination in both developing and diseased brains. T2-, and T2*-weighted multi-echo data have been proposed to estimate MWF in the human brain. Even for the simple two pool signal models of myelin and non-myelin associated water, number of dimensions of the parameter space for obtaining MWF estimates remains high, making parameter estimation challenging. The aim of this research is to improve the accuracy of brain myelin mapping using a novel joint inversion imaging method that combines data from multiple contrasts in a single optimization process.

Introduction

Brain myelin content is a significant imaging diagnostic biomarker that can be measured with MRI by estimating the amount of myelin water (MW), i.e., the amount of water trapped in the myelin sheath. MW has shorter relaxation times than the surrounding free water (FW). The myelin water fraction (MWF = amount of MW divided by the total amount of water) is derived from the peak areas of relaxation spectra, assuming an upper threshold for the MW component. Commonly multi-exponential relaxometry methods are employed, using either multi-echo spin-echo (ME-SE) or multi-echo gradient-echo (ME-GE) experiments^1,2. The underlying inversion problem is ill-conditioned; it is susceptible to noise and measurement inadequacies, potentially leading to less accurate and more biased estimated parameters³. There are numerous methods available for extracting multi-exponential parameters from a time series. Non-negative least squares (NNLS) with regularization is the most commonly used method for MWF estimation^4,5. Here we introduce a Levenberg-Marquardt (LM) inversion algorithm to jointly recover T2- and T2*-distribution, which provides MWF maps with fast computational speed and stable solutions. The method is inspired from geophysical joint inversion, where several ill-posed inverse problems are frequently jointly solved to provide stabilized and more accurate results^6,7. In this research, we use the multi-exponential signal model and extended phase graph (EPG) formalism⁸ to generate synthetic ME-GE and ME-SE signal decay curves, respectively, with realistic myelin water composition. MWF analysis was performed with NNLS and the new LM approach, and results were compared with the ground truth. We address the effects of noise on the computation in scenarios with low to high MWF.

Materials and Methods

To simulate decay curves with continuous distribution comparable to white matter, two Gaussian distributed relaxation compartments were used as shown in Figure 1 (Input Model). For myelin imaging, the distribution are thought to have two peaks, a smaller peak for MW and a larger peak representing FW. Using the EPG formalism and a simple exponential signal decay model, synthetic ME-SE and ME-GE signals were derived from the T2 & T2* distribution, respectively (cf. Figure 1, Synthetic Signals). We normalize (set integral to one) the distributions so that the integral of the fast component is directly the MWF. We can then parametrize the distributions with 9 parameters: means and standard deviations of both peaks for both distributions, and the integral of the fast component (same value for both distributions). We propose a method to estimate these parameters by forming a non-linear minimization problem (Figure 1, Inversion) with an objective function based on the difference between the "observed" (synthetic) and predicted decay curves. We iteratively solve the non-linear problem with a local-search optimization method, the Levenberg-Marquardt algorithm. At each iteration, we first perform the computation of both decay curves as described before, and then compute the Fréchet derivatives. We solve the linearized problem as a classical linear least-squares problem with bounds on the variables. We added a damping term to the cost function for more stability. In order to simulate complete imaging data, we synthesized ME-GE and ME-SE data for a given MWF distribution in the brain, which was taken from a published MWF atlas⁹. We add Rician distributed noise to the simulated data and perform voxel-vise MWF analysis with common NNLS and the newly proposed LM approach for individual inversion of ME-SE and ME-GE data, respectively. Finally, the LM joint inversion is performed on the combined data. The inversion is performed independently on each voxel and outputs the perturbations of the model parameters to be added to the initial model. From the final parameters obtained after the inversion, the last one gives us directly the MWF. We compute the deviation to the input MWF (ground truth) and its Root-Mean-Square-Error (RMSE) across all voxels, which serves as a metric for the accuracy of the inversion. An Apple MacBook with M1 Pro Silicon chip, 32GB RAM, and with Python 3.8 installed is used for implementing the whole pipeline. For the inversion of the MWF map, the LM approach (less than 120 seconds) was computationally less expensive as compared to NNLS (around 300 seconds).

Results

Figure 2 presents the results of the estimated MWF maps using NNLS and LM for individual inversion of GE and SE data, both for SNR values of 20 and 50. As expected, RMSE decreases when SNR increases. The LM based method generally shows better performance (i.e lower RMSE) than NNLS. The MWF maps obtained through joint inversion (Figure 3) further reduce the RMSE as compared to individual inversion with LM. When using GE data with SNR=20 and SE data with SNR=50 (Fig. 3, bottom row) the joint inversion result is nearly as good as for the case SNR=50 in both datasets.

Discussion

When compared to the commonly used regularized NNLS method, the LM approach results in reduction of absolute error for the estimation of MWF using either ME-GE or ME-SE data. The LM joint inversion of the combined data strongly reduces the error compared to the individual inversion. These results were obtained on synthetic data with realistic SNR values, which enable a direct comparison to the ground truth MWF. The next step is to test the approach on real MR data.

Acknowledgements

This work received financial support from the Helmholtz Association, Initiative and Networking Fund, funding code ZT-I-PF-4-006 (Helmholtz Imaging Project "JIMM”).