Yudu Li^{1,2}, Rong Guo^{1,2}, Yibo Zhao^{1,2}, Yang Chen^{1,3}, Bryan Clifford^{1,2}, Tianyao Wang^{4}, Chenyan Wang^{1,5}, Yiping Du^{6}, Yao Li^{6,7}, and Zhi-Pei Liang^{1,2}

Quantitative mapping of myelin water fractions (MWF) can substantially improve our understanding of the progression of several demyelination white matter diseases such as multiple sclerosis. While MWF can be determined from both T2-weighted and T2*-weighted imaging data, it is much faster to collect T2*-weighted imaging data. However, estimation of MWF from T2*-weighted imaging data using a multi-exponential component model is an ill-conditioned problem whose solutions are often very sensitive to noise and modeling errors. In this work, we address this problem using a new model-based method. This method is characterized by: a) absorbing the spectral priors using the Bayesian-based statistical framework, and b) absorbing the spatial priors using a reproducing kernel based model. Both simulation and experimental results show the proposed method significantly outperforms the conventional parameter estimation methods currently used for MWF estimation.

A common model used for MWF estimation decomposes
the T2*-weighted
imaging data into three exponential components^{3-5}:

$$\hspace{18em}\rho(\boldsymbol{x}_p,t)=\sum_{l=1}^3A_{\ell}(\boldsymbol{x}_p)e^{-t/T_{2,\ell}^{*}(\boldsymbol{x}_p)}e^{i2{\pi}f_{\ell}(\boldsymbol{x}_p)t},\hspace{18em}\quad\quad(1)$$

which correspond to the water within the myelin sheath, myelinated axons and all others tissues. $$$A_{\ell}(\boldsymbol{x}_p)$$$, $$$T_{2,\ell}^{*}(\boldsymbol{x}_p)$$$ and $$$f_{\ell}(\boldsymbol{x}_p)$$$ denote the amplitude (related to concentration), $$$T_2^*$$$ relaxation time constant and frequency shift of each component at location $$$\boldsymbol{x}_p$$$ ($$$p=1,2,...,P$$$). Accurate signal decomposition based on Eq. (1) is challenging due to ill-conditonedness. Here, we address this issue by incorporating both spectral and spatial prior information.

**Incorporation of spectral priors**

Spectral priors are incorporated using the Bayesian-based statistical framework. More specifically, we estimate the nonlinear spectral parameters (i.e., $$$\{T_{2,\ell}^{*}\}$$$ and $$$\{f_{\ell}\}$$$) in two steps.

Firstly, we obtain the initial estimates of $$$\{T_{2,\ell}^{*}\}$$$ and $$$\{f_{\ell}\}$$$ using the conventional parameter estimation methods (we denote the solution as $$$\{\bar{T}_{2,\ell}^{*}\}$$$ and $$$\{\bar{f}_{\ell}\}$$$ ) and then approximate their empirical distributions using a Gaussian model:

$$\text{Prob}(T_{2,\ell}^*){\approx}\mathcal{N}(\mu_{T_{2,\ell}^{*}},\sigma_{T_{2,\ell}^{*}}^{2})\quad\text{and}\quad\text{Prob}(f_{\ell}){\approx}\mathcal{N}(\mu_{f_{\ell}},\sigma_{f_{\ell}}^{2}).$$

Secondly, we re-estimate the nonlinear parameters, incorporating the empirical prior distributions $$$\text{Prob}(T_{2,\ell}^*)$$$ and $$$\text{Prob}(f_{\ell})$$$ via the Bayesian-based statistical framework:

$$\tilde{T}_{2,\ell}^{*}=\arg\min_{T_{2,\ell}^*}\sum_{p=1}^{P}||\bar{T}_{2,\ell}^*(\boldsymbol{x}_p)-T_{2,\ell}^*(\boldsymbol{x}_p)||_2^2+\lambda_{T_{2,\ell}^*}||T_{2,\ell}^*(x_p)-\mu_{T_{2,\ell}^*}||_2^2$$

$$\tilde{f}_{\ell}=\arg\min_{f_{\ell}}\sum_{p=1}^{P}||\bar{f}_{\ell}(\boldsymbol{x}_p)-f_{\ell}(\boldsymbol{x}_p)||_2^2+\lambda_{f_{\ell}}\sum_{p=1}^P||f_{\ell}(x_p)-\mu_{f_{\ell}}||_2^2.$$

**Incorporation of spatial priors**

Once the nonlinear parameters have been determined, we
estimate the linear coefficients $$$\{A_{\ell}(\boldsymbol{x}_p)\}$$$ for all the voxels jointly, incorporating
spatial priors. To this end, we use a kernel-based model to represent $$$A_{\ell}(\boldsymbol{x}_p)$$$, leveraging our previous
work in this area^{6}. More specifically, we assume $$$A_{\ell}(\boldsymbol{x}_p)$$$ can be viewed as a nonlinear function of a set
of low-dimensional features $$$\boldsymbol{z}_p{\in}{\mathbb{R}^n}$$$ derived from prior images (e.g., T1W images). While
such a nonlinear function is often highly complex in practice, we can linearize
it in a high-dimensional transformed space using the “kernel trick” in machine
learning:

$$\hspace{22.3em}A_{\ell}(\boldsymbol{x}_p)=\sum_{q=1}^P\alpha_{\ell,q}k(\boldsymbol{z}_q,\boldsymbol{z}_p),\hspace{22.3em}\quad(2)$$

where $$$k(\cdot,\cdot)$$$ is some kernel function. Substituting Eq. (2) into Eq. (1) leads to the kernel-based representation for the T2*-weighted imaging data:

$$\hspace{16.5em}\rho(\boldsymbol{x}_p,t)=\sum_{\ell=1}^3\left\{\sum_{q=1}^P\alpha_{\ell,q}k(\boldsymbol{z}_q,\boldsymbol{z}_p)\right\}e^{-t/T_{2,\ell}^{*}(\boldsymbol{x}_p)}e^{i2{\pi}f_{\ell}(\boldsymbol{x}_p)t}.\hspace{16.5em}\quad(3)$$

To determine the $$$\{\alpha_{\ell,q}\}$$$ in Eq. (3), we choose the radial Gaussian
kernel function and fit the kernel-based model to the measurements by solving a
linear least-squares optimization problem that takes into account all the
voxels jointly. Once the $$$\{\alpha_{\ell,q}\}$$$ have been estimated, we can reconstruct $$$A_{\ell}(\boldsymbol{x}_p)$$$ based on Eq. (2), which will yield the desired
MWF maps^{3-5}.

The performance of the
proposed method has been evaluated using both simulated and experimental data.
In our Monte-Carlo simulation study, the simulated data were generated based on
Eq. (1) using published literature values^{3} for the spectral
parameters and in vivo levels of additive Gaussian white noise. Figure 1 shows the estimation standard deviation from the Monte-Carlo study and the MWF maps estimated from
one of the 40 realizations. As can be seen, the MWF estimates obtained using
the proposed method show significantly reduced errors and estimation variances.
In our experimental study, in vivo data were collected on a 3T scanner
using a multi-gradient-echo sequence. Two sets of representative results are
shown in Fig. 2 and 3, one from a healthy subject (Fig. 2) and the other from a
patient with a chronic stroke (Fig. 3). As can be seen, the MWF maps produced
by the proposed method show significantly reduced artifacts and estimation
variations over those obtained using the conventional scheme, which is consistent
with those from our simulation study.

[1] MacKay AL, Laule C. Magnetic resonance of myelin water: an in vivo marker for myelin. Neural Plast. 2016;2(1): 71-91.

[2] Alonso‐Ortiz E, Levesque IR, Pike GB. MRI‐based myelin water imaging: A technical review. Magn Reson Med. 2015;73(1): 70-81.

[3] Du YP, Chu R, Hwang D, et al. Fast multislice mapping of the myelin water fraction using multicompartment analysis of T2* decay at 3T: A preliminary postmortem study. Magn Reson Med. 2007;58(5):865-870.

[4] Hwang D, Kim DH, Du YP. In vivo multi-slice mapping of myelin water content using T2* decay. NeuroImage. 2010;52(1):198-204.

[5] Nam Y, Lee J, Hwang D, et al. Improved estimation of myelin water fraction using complex model fitting. NeuroImage. 2015;116:214-221.

[6] Li Y, Liang ZP. Constrained image reconstruction using a kernel+sparse model. Proc Intl Soc Magn Reson Med. 2018;p. 657.

Figure 1: Simulation results from a Monte-Carlo
study: (a) MWF maps estimated in one of the 40 realizations, and (b) estimation standard deviation (STD) maps (within the gray matter and white matter only)
normalized by the true MWFs. The left column shows the results obtained by the
conventional parameter estimation method and the right columns shows the results produced by the
proposed method. The averaged STD for both methods
are also shown in red text.

Figure 2. In vivo experimental results of the
MWF maps from a healthy subject at three representative slices across the brain,
estimated using: (a) the conventional parameter estimation method, and (b) the
proposed method. Note that the MWF maps from the proposed method show
significantly reduced estimation variations and much better definition of tissue
boundaries, which is consistent with the results from our simulation study.

Figure 3. In vivo experimental results from a
stroke patient at two representative slices containing a chronic stroke,
including (a) structural images, (b) MWF maps obtained by the conventional parameter
estimation method, and (c) MWF maps produced by the proposed method. The
location of the chronic stroke is indicated by red arrows. As can be seen, the
MWF maps estimated by the proposed method show a significant reduction in MWF
within the stroke, which is barely observable in the results obtained from
the conventional method.