Improved magnetic resonance myelin water imaging using multi-channel denoising convolutional neural networks (MCDnCNN)

Introduction

Myelin water fraction (MWF) derived from quantitative myelin water imaging (MWI), is a powerful and important magnetic resonance imaging (MRI) method for studying and diagnosing neurological and psychiatric diseases such as multiple sclerosis (1), schizophrenia (2), Alzheimer disease (3) and stroke (4,5). MWF maps are commonly calculated based on T₂/T₂* imaging data acquired using a Carr-Purcell-Meiboom-Gill (CPMG) sequence or its accelerated version, i.e., the gradient and spin echo (GRASE) sequence, in which the T₂ value of the multicomponent water is analyzed (6). Recently, a multiple gradient recalled echo (mGRE) sequence (7-9) was developed as an acquisition method for MWI data, in which T₂* decay signals are analyzed for the entire brain. Compared with the CPMG or GRASE sequence, the mGRE sequence has the advantages of faster acquisition time, lower specific absorption rate, shortened first echo time (TE) and echo spacing.

MWF can be measured from mGRE through a complex three-pool model (3CCT₂*) with each pool representing a different water component and a frequency offset term added to each of the three components to minimize extra frequency shifts that may possibly be introduced. Such frequency shifts depend on the angle of the underlying neural fibers orientation relative to the main magnetic field B0 (9-12). However, the sensitivity and instability of multiexponential decay analysis to noise (13), a well-known difficulty in MWF estimation (14-16) that will lead to inaccuracies in deriving MWF estimates, should be taken into consideration. In particular, the sensitivity also makes it more difficult to estimate MWF mapping at higher spatial resolutions, although higher resolutions may help to treat partial volume effect, detect small demyelization lesions and practice early diagnosis. The situation can be eased if the quantitative susceptibility mapping (QSM) pipeline is employed, for overcoming non-local phase contrast and orientation dependence (12). Then the challenges would turn to be the issue of signal-to-noise ratio (SNR), due to the reduced size of the voxels. To maintain necessary and adequate accuracy and precision, noise reduction filters should be applied in the stage of data postprocessing (8,16,17). However, the conventional linear filtering has a number of disadvantages, including introduction of possible bias that may lead to inaccurate parameter estimation and obscuration of small-scale structures due to partial volume effect (18). Although some advanced filters have been applied to improve MWF mapping, the performance of these advanced filters will decline when dealing with data containing moderate levels of SNR, as shown by their results (8,10,12).

A new denoising method was recently introduced. It was a feed-forward multi-channel denoising convolutional neural networks (MCDnCNN), which was originally proposed for removing noise from 3D volumes of MRI data (19). In contrast to a direct estimation of noise-free truth image (i.e., the latent clean image), MCDnCNN adopts a residual learning strategy to isolate noise rather than the latent clean image from its noisy observation (20). In other words, MCDnCNN has been designed so to predict the noise at each voxel of the raw image. The MCDnCNN has a training mechanism that is faster, more stable and more robust for denoising, compared to using the traditional 2D denoising network (19). Up to date, MCDnCNN has outperformed all other tested methods in study, including non-local means and discrete cosine transform, in terms of noise reduction and detail preservation (19). Moreover, MCDnCNN also performs much faster than other state-of-the-art methods for noise reduction (19,20).

Our hypothesis was that MWF calculated from mGRE data using the MCDnCNN models corresponding to the specific noise level (denoted as appropriate MCDnCNN models) should be both qualitatively and quantitatively improved over that calculated from noisy MRI data or inappropriately MCDnCNN-denoised mGRE data, as the appropriate MCDnCNN models should be able to estimate the noise more accurately. Our procedure first trained the five different MCDnCNN models corresponding to four specific noise levels and a mixed noise levels using real-world in vivo brain imaging data acquired at the first TE of the mGRE sequence. We also simulated data in the brain by applying the hollow cylinder fiber model (HCFM), recently proposed by Wharton et al. (21), to generate complex signals of the magnitude and phase of the multicomponent mGRE data at five different noise levels. Next, we applied the five trained MCDnCNN models to denoise the simulated signals and in vivo data. Lastly, we analyzed these denoised signals of both the simulated and in vivo data using the 3CCT₂* fitting models (9,12). We examined the results against published in vivo MWF results to assess the performance of our proposed method.

Methods

The main technique employed for evaluating the improvement in estimating MWF in this study was MCDnCNN, with the initial susceptibility in the data estimated using a method called sparse linear equation and least-squares (iLSQR) (22) based on the 3CCT₂* model. To highlight how MCDnCNN may improve the accuracy of estimating MWF maps, we adopted two procedures for obtaining MWF maps, one using MCDnCNN in the procedure and the other not using it. In the former procedure, MCDnCNN was applied to obtain denoised magnitude images from original magnitude images. Meanwhile, iLSQR was used to obtain QSM (i.e., the magnetic susceptibility maps) from the original phase images. The denoised images and the QSM were synthetized to estimate the decay curves and then to evaluate MWF maps using 3CCT₂*. The latter procedure did the same except for not including a step of using MCDnCNN (Figure 1). We will elaborate the details in the procedure of using MCDnCNN (Figure 2).

Figure 1 The general processing flow evaluating MWF in our protocol. The original magnitude images were first denoised using MCDnCNN; and its corresponding phase images were used to calculated QSM maps (magnetic susceptibility maps). The denoised images and QSM maps were used together to estimate the T₂* decay curves and finally to compute the MWF maps based on the 3CCT₂* model. QSM maps were useful in overcoming non-local phase contrast and orientation dependence. The flow denoted by the dark and red arrows represents the procedure of evaluating MWF using MCDnCNN (the first procedure). The flow denoted by blue and red arrows represents the procedure of evaluating MWF without using MCDnCNN (the second procedure). MCDnCNN, multichannel version of the denoising convolutional neural networks; iLSQR, initial susceptibility estimated by sparse linear equation and least-squares; QSM, quantitative susceptibility mapping; 3CCT₂*, complex three-pool model; MWF, myelin water fraction.

Figure 2 The network structure of the ten layers of the MCDnCNN learning model. The first layer is the input layer (purple) which consists of a convolution layer (blue). Each of the eight hidden layers (yellow) consists of a BN layer (green) and a convolution layer (blue). All convolutional layers are to be activated by ReLU (red). The last layer is the output layer (blue) which consists of a convolution layer (blue). The lower part in this figure illustrates the internal structures of the input and hidden layers, respectively. The note of ‘3×3×64, 64’ on the hidden layers denotes that each hidden layer processes 64 kernels each with an image window of dimensions =3×3×64. S−2 to S+2, a 5-slice block from slice location −2 to +2 centering the current slice S (the first to the fifth slice of the five-slice block); 64, number of features; 3×3×64/3×3×5, size of convolution kernel; Conv, convolution layer; ReLU, rectified linear unit; S, the third (middle) slice of the five-slice block; MCDnCNN, multichannel version of the denoising convolutional neural networks; BN, batch normalization.

Data acquisition

The T₂* imaging data and T1 weighted images were acquired on a 3-T MRI scanner (Siemens, Erlangen, Germany) with 20-channel radio frequency head-coil at Shanghai Key Laboratory of Magnetic Resonance. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). The study was approved by the local Ethics Review Board of East China Normal University (IRB#HR227-2019). Informed consents were signed by 26 healthy adult volunteers (mean age 22.81±2.47 years). The imaging sequence was a 24-echo 3D gradient recalled echo scan with bipolar echo readout using the following acquisition parameters in axial plane: first echo time (TE1) =2.6 ms and echo spacing time (ESP) =1.50 ms, repetition time (TR) =62 ms, flip angle =15°, slice thickness =2.0 mm, voxel size =0.9×0.9×2.0 mm³, matrix size =288×288×64, bandwidth =870 Hz/voxel and phase encoding direction acceleration (GRAPPA) =2, field of view (FOV) =256×256×128 mm³. Duration of data acquisition was 10 min 42 s. The dimensions of the final images were 288×288×64 without zero-padding. T1-weighted images were acquired using a magnetization-prepared rapid gradient echo sequence (MPRAGE). The MPRAGE parameters were: 192 sagittal slices, TR/TE/inversion time (TI) =2,530/2.98/1,100 ms, flip angle =7°, slice thickness/gap =1/0.5 mm, voxel size =1.0×1.0×1.0 mm³, matrix size =256×256×192, FOV =256×256×288 mm³. Duration of data acquisition was 6 min 03 s. Foam pads were placed around each subject’s head during each scan to minimize head motion. When subject motion was observed, the scan was repeated.

MCDnCNN

The MCDnCNN proposed by Jiang et al. (19) was developed based on a residual network for denoising 3D images. The goal of MCDnCNN was to compute the noise-free image of a noisy input image. The proposed MCDnCNN is a complicated multi-layer residual learning network (23), which explicitly fits a residual mapping using its deeper network containing multiple layers and is easier to optimize all the weights than conventional neuron learning (Figure 2). The structure of MCDnCNN consists of one input layer of convolution with a rectified linear unit (ReLU) (24) as an activation function (19,20), eight layers of convolution with batch normalization (BN) (25) (a normalization function for each mini-batch of the training to speed up training), and finally one output layer of convolution. Residual learning and BN are introduced to facilitate network optimization and to speed up the training process (23,25). The network implicitly removes the latent clean image through the 8 hidden layers so that an image of noise can be estimated, and the noise-free image can be generated by subtracting the estimated noise image from the original image.

This network (Figure 2) does not have any max-pooling layer, which is a feature reduction facility for identifying the most dominant features (26). Therefore, the output layer would generate an image of the same size as that of the input. For the input layer of convolution with an ReLU, we used 64 different kernels of size 3×3×5, whose weights were assigned with randomly generated real-values at the start of training. The input volume was split into several 3D blocks of size X × Y × 5 with a stride of 25 in the XY-plane and 1 in the Z-direction, and each individual block was used as input to the network to take into consideration the neighborhood context of the current slices. Thus, the first layer would result in 64 images of the same size for each of the X × Y × 5 image patch. Each of the 8 hidden layers employed a similar set of 64 individual 3×3×64 kernels while the last output layer employed a single 3×3×64 kernel to continue processing the data resulted from the previous layer. Residual learning formulation and BN, as two key features (20) of this network, were incorporated to speed up training as well as to boost the denoising performance. During the training process, the averaged mean square error between the desired residual image (the noise that we superimposed) and the estimated ones (output from the network) from the noisy input was applied as the loss function of the networks to update and optimize the parameters of the MCDnCNN. The original image before noise was superimposed or the original noise free volume should be restored after the noise at every slice in the entire volume was estimated.

We trained the MCDnCNN model at four specific noise levels and a mixed noise level. The noise took use of the standard deviation (SD) of the background in the image acquired at the last TE while the signal took use of the mean intensity in homogeneous regions of white matter in the image acquired at the TE1. The homogeneous regions were manually selected within one white matter structure visually identified with the biggest size in the image. The noise level was then computed as a ratio of noise divided by signal, basically a reverse of SNR for the entire mGRE dataset. We assumed in this scenario that TE1 images and the last TE images of the mGRE data were clean images. Then we added a specific level of Rician noise to both the first TE and last TE images to get the noised first TE and noised last TE images, and next evaluated the noise level based on the two noised images. The Rician distribution is a function of the underlying magnitude. Thus, adding Rician noise is defined as:

R=(r+δ×s)2+(δ×s)2;s=randn(0,1)[1]

where δ is defined as the ratio between the noise value to be added and the maximum intensity of the image (such that δ =1% means the noise value 0.01 if the maximum intensity of image is 1), r and R represent the images (intensity-normalized) before and after adding the noise, respectively. “randn(0,1)” is a function that generates random values following a normal distribution (mean =0, SD =1). For instance, we added 1% (δ =1%) Rician noise to both the first and last TE images to get the noised TE1 images and noised last TE images, and thus obtained Level1 defined by the ratio of the noise (in the image acquired at the last TE) divided by the mean signal intensity (of the image acquired at the first TE). The mGRE images at TE1 from human volunteers, which was added with the noise at the specific noise level, were used for training the MCDnCNN model. The SNR of the images at TE1 is typically highest before noise is added to the images, thus TE1 images were chosen for training the networks while images acquired at other echoes were not used for the training. Images of the 15 out of the 26 participants were randomly selected to add noise at a specific noise level (e.g., Level1, Level2, Level3, Level4, or LevelMix) for training the neural network (10 for training and 5 for validation). Noise was not added to the rest 11 participants, and they were used as real-world datasets to evaluate the models. In the following text, Level1 (noise level =0.0065), Level2 (noise level =0.0113), Level3 (noise level =0.0205), Level4 (noise level =0.0293) and LevelMix denote the first and last TE images with added Rician noise (δ) at levels of 1%, 3%, 5%, 7% and also 1–7% (mixed), respectively; and Delevel1, Delevel2, Delevel3, Delevel4 and DelevelMix denote the MCDnCNN denoising models corresponding to the noise levels, respectively. It should be noted that LevelMix actually added noise to each image with noise levels ranging from 1% to 7%. And we randomly picked the same number of patches with the other four noise levels to train the DelevelMix (19).

For training the models to denoise Rician noise at a noise level, patch size X × Y × 5 was set to 60×60×5, and thus for the 15 training datasets with the sliding window, about 150,000 patches were obtained to train the corresponding model. The Adam’s gradient-based stochastic optimization algorithm was used to update the weights (27), which took the averaged mean square error as our cost function. To speed-up the optimization process, the windowed sub-images were fed to the training model using a so-called mini-batch strategy, a well-established method for addressing stochastic optimization problems (28). The mini-batch size was set to 64 patches as allowed by our hardware to take full use of the GPU. The MatConvNet software package (29) was employed for the training of the MCDnCNN models. The training processes converged rapidly within only about 50 iterative epochs, in which the learning curves decayed exponentially and quickly moved from along 1e⁻¹ to 1e⁻⁴. Each training took about 8 h on a personal computer equipped with an NVIDIA GeForce GTX 1050 when the computer was fully dedicated. Training to meet clinical needs on the fly was not practical.

MWF mapping

Signal processing and QSM

After applying channel-wise inverse Fourier transform to obtain the coil magnitude and phase images, followed by combining coil images to obtain robust echo magnitude and phase images, the mGRE volumes at odd-numbered echoes of TE1, 5, 9, 13, 17 and 21 in the sequence with an interval of four ESP’s (i.e., 1.5 ms × 4=6 ms) were used for QSM calculation. The ESP was chosen so just to follow popular acquisition protocols for QSM (12). The magnitude image of the first echo was adopted to form a mask using the BET tool provided as part of FSL software package (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/). The mask was then binarized. STI Suite V3.0 (https://people.duke.edu/_cl160/) was used to calculate QSM and by the way handled susceptibility maps. In particular, 3D V_SHARP (for removing the background phase) and iLSQR (22) in the STI Suite were used to generate 3D susceptibility maps.

Signal synthetization

To reduce the influence of non-local magnetic fields, we computed quantitative susceptibility maps based on the phase information of the tissue. The susceptibility maps were subsequently used to estimate frequency shift to overcome drawbacks known in the procedure working directly with tissue phase (12) using:

Δf=χγB02π[2]

where χ denotes the magnetic susceptibility, whose value can be estimated using the STI Suite pipeline, γ is the gyromagnetic ratio of hydrogen and B0 =2.89 T (for the nominally 3 T Siemens scanner) is the strength of the static magnetic field (30-32). Eq. [2] was particularly tailored for the current study as it may minimize the impact from the background field and preserve the local field introduced by the local underlying tissue and it was regarded as a denoising operation for the phase (12,33). Subsequently, the voxel-wise corrected version of the complex signal S_m(t) at each TE time point t can be obtained using magnitude A_m and frequency shift Δf_m:

Sm(t)=Am(t)e-j2πΔfmt[3]

where, m denotes the ordinal number of TE, i.e., TE1, TE2, …, and so on and forth. In accordance with previous studies (12), compartments of the complex signal were modelled using the complex three-pool model {Eq. [4]} for the mGRE T₂* data (3CCT₂*):

S(t)=[AMWe−(1/T2,MW∗+j2πfMW)t+AAWe−(1/T2,AW∗+j2πfAW)t+AEWe−(1/T2,EW∗+j2πfEW)t]ejφ0[4]

where MW = myelin water, AW = axonal water, EW = extracellular water; A’s and f’s respectively represent signal magnitudes and frequency offsets in each pool, and φ₀ is the initial phase offset. Following the outcomes of our experiments and those from previous published reports, the searching boundaries and initial values of each variable in Eq. [4] were set according to those in the literature (9,12) assuming no background components in the frequency terms. The searching boundaries of A_MW, A_AW and A_EW were from 0 to 2 times of the magnitude of the complex signal at the first echo (|S1|) while the initial values were set to 0.1×|S1|,0.3×|S1| and 0.6×|S1|, respectively. The range of T₂* of MW was limited to between 3 and 25 ms while the initial value was set to 10 ms. The ceiling and floor values of T₂* of axonal water and extracellular water were 25 and 150 ms, while the initial values were set to 64 and 48 ms, respectively. The range of f_MW was −75 to +75 Hz while the initial value was set to 5 Hz. The searching boundaries of f_AW and f_EW were −25 to +25 Hz and the initial values were set to 0. The upper and lower values of the last term, φ₀, were −π to +π while the initial value was the angle of the signal of the first echo. All the searching boundaries and initials were consistent with those in the reference (12).

After the magnitudes of each myelin component, i.e., A_MW, A_AW, A_EW, in Eq. [5] were estimated, the MWF can be calculated as the fraction of MW over the total water signal amplitude:

MWF=AMW/(AMW+AAW+AEW)[5]

The MWF values were quantified based on imaging data before and after denoising, respectively.

Data simulation and corresponding numerical analysis

A HCFM was used in the simulation. T₂* components were generated for MW , axonal water, and extracellular water, with relaxation times for data simulated at 3 T equal to 10, 64, and 48 ms (9), respectively. The percentages of the 3 components in the brain were set as 12%, 38% and 50% (9,34). A total of 290,000 MWF values of these components, ranging randomly in the scope of 5–18% (MWF =5–18%), were simulated to reflect a range of myelin fraction from low to high. The values of the three pools, MW, axonl water and extracellular water varied with MWF and they summed up to 100%, mimicking the real cases in the in vivo data. We calculated a white matter mask of MWF map using an in vivo brain dataset and assigned the simulated MWF values to the voxels within the masked area. We sorted and relocated the simulated MWF values within the masked area so that the distribution of the MWF values was consistent with that of the in vivo brain. The manufactured MWF map was to be used as a ground truth later in our procedures for validating the calculated MWF maps generated based on MCDnCNN-denoised data. The frequency offsets for each of the three-pool components were calculated using a simplified 1D expression of the HCFM model (21) with an orientation setting of five different angles ranging from 0 to π/2 with an incremental step of π/8.

The HCFM was used to predict average frequency perturbations (Δf) in each radial compartment [see (21) for a detailed derivation]. According to HCFM, the frequency shift of myelin water (f_MW) can be expressed as follows:

fMW=γB02π(χI2(cos2θ−13)+E+χA2(−13+sin2θ(14−32ri2ro2−ri2ln(rori))))[6]

where γ represents the gyromagnetic ratio and θ the angle between the orientations of the local fiber bundles and the magnetic field B0. The following parameters were used to calculate the frequency offset of MW: isotropic susceptibility (χ_I) =−100 ppb; anisotropic susceptibility (χ_A) =−100 ppb; chemical exchange (E) =20 ppb; g-ratio (r_o/r_i where r_o: outer radius and r_i: inner radius) =0.8; and B0 =3 T (9,35). Complex signals were generated for each orientation with the 24 echoes, TE1 =2.6 ms and ESP =1.5 ms. Rician noise (δ) was added only in magnitude images at 5 different levels 1%, 3%, 5%, 7% and also 1–7% (mixed), which corresponded to Level1 (noise level =0.0065), Level2 (noise level =0.0113), Level3 (noise level =0.0205), Level4 (noise level =0.0293) and LevelMix, respectively, as mentioned earlier.

The simulated signals were then processed and analyzed following the procedures as previously described in the text, resulting in two versions, i.e., denoised by MCDnCNN or not denoised. Simulated data at every noise level were denoised by all five MCDnCNN models. For instance, the data at noise Level1 were denoised by Delevel1, Delevel2, Delevel3, Delevel4 and DelevelMix, respectively. To determine the acceptable level of mismatched noise, we also conducted an additional experiment. We chose to use the simulation data at Level2 (noise level =0.0113) and Level3 (noise level =0.0205) with an orientation angle of π/2 (the myelinated fibers were perpendicular to B0 direction). The increments of adding Rician noise were ±0.25%, ±0.5% and ±0.75%, and the corresponding noise levels were also calculated, respectively. Delevel2 and Delevel3 were used to denoise the corresponding noisy images, respectively. All the results were then fitted using the complex three-pool model to calculate different versions of MWF maps. To demonstrate the effect of using MCDnCNN, the mean and SD of the absolute error (AE) between the calculated MWF (either before or after the denoising) and the ground truth MWF were computed. All simulations were performed within MATLAB (Matlab R2018b, The MathWorks Inc., Natick, MA, USA).

Analysis of in vivo data

T₂* decay curves were first corrected using the corrected phase information and the original data. MWF maps were then calculated using the complex three-pool model. The same also applied to the mGRE data that were already denoised using MCDnCNN, so that another set of T₂* decay curves and MWF maps were obtained. The 3CCT₂* model is sensitive to noise. To explore the sensitivity of 3CCT₂* to noise and the effectiveness of the appropriate MCDnCNN models to obtain the MWF values, one original in vivo dataset was picked randomly. In the process of denoising, the noise level of the data was calculated and the five MCDnCNN models were applied to this dataset. At last, the four denoised T₂* decay curves were obtained and compared with the T₂* decay curves estimated based on the data without being denoised. To visualize the results in a more intuitive and informative way, we randomly sampled thirty voxels for reference from the 3D mGRE maps to illustrate the sensitivity and effectiveness. Meanwhile five versions of the denoised MWF were computed by 3CCT₂* and compared with the MWF values calculated based on the data without being denoised (30 voxels for reference, two sample t-test, P<0.05, uncorrected).

All 11 subjects were repeated following the same procedures that calculated noise level and computed MWF values before and after denoising. Also, to confirm that the appropriate MCDnCNN models did not affect image contrast, we also analyzed contrast-to-noise ratio (CNR) of the data before and after denoising (36). Five regions-of-interest (ROIs) in the brain of each of the 11 subjects were selected to perform statistical comparisons of the MWF values in terms of mean and SD (two sample t-test, P<0.05, uncorrected). To examine the effectiveness of our method, the results were quantitatively validated against those published in the literatures (10,12,37).

Results

Simulation experiments

MWF maps derived from the images denoised by appropriate MCDnCNN models [e.g., Level3 (δ =5%) should be denoised by Delevel3] exhibited that edges and small structures were well preserved and using appropriate MCDnCNN models may recover most of the structural details from noisy images (Figure 3). Checking against the ground truth MWF map, we may see that some of the structures were not completely recovered in the maps that were calculated based on the data without being denoised (e.g., the areas in the green and red rectangles in Figure 3). However, the structure in the Level3 data denoised by Delevel3 was appropriately preserved or restored. Nevertheless, if denoised by inappropriate MCDnCNN models, for example, by Delevel1, Delevel2, Delvel4 or DelevelMix, the structures in the ground truth indicated by the green and red rectangles were hardly recovered and the corresponding MWF values appeared to be overestimated, too. Furthermore, the mean and SD of AE values derived from images denoised using the appropriate MCDnCNN models were reduced, compared with those derived from the original data before denoising (Tables 1,2). In particular, taking the Level3 data (orientation angle =0) as an example, we may see that the mean and SD of AE values based on the data denoised by Delevel3 (an appropriate model) were 0.09739 and 0.07907 whereas those derived directly based on the original data were higher at 0.1177 and 0.09360, respectively. In contrast, the mean and SD of AE values derived from images that were denoised using inappropriate MCDnCNN models exhibited systematic overestimates and overall larger SD values (Figure 3, Tables 1,2). For example, when the Level2 data (orientation angle = π/8) were denoised inappropriately using Delevel1, the mean and SD of AE values derived from Delevel1 were 0.08911 and 0.07253 whereas those derived from its appropriate Delevel2 model were 0.07346 and 0.06010, respectively. Again, it showed that using inappropriate models would inappropriately overestimate MWF values as the corresponding mean values of AE were obviously high.

Figure 3 A comparison of the MWF maps based on the simulation data of Level3 denoised by the five models. A typical slice at the same slice location of the five results is displayed. The top row is the ground truth of the MWF map. Figures in panel A from the second row to the last row (including the two-small zoomed-in insets) are the MWF maps at the same slice location of the 6 results (2nd row: derived from the original raw data without denoising; 3rd–7th rows: derived from the data denoised by the 5 models of Delevel1, Delevel2, Delevel3, Delevel4 and DelevelMix, respectively). Figures in panel B from the second row to the last row are AE maps of MWF at the same slice location corresponding to those in panel A. Figures in panel C from the second row to the last row are the corresponding AE histograms, with the mean and SD of AE values noted in red. Colored bars represent the MWF values. Delevel, the denoising model corresponding to the noise level; MWF, myelin water fraction; AE, absolute error; SD, standard deviation.

In addition, the mean and SD of AE values derived from data denoised using the general MCDnCNN model DelevelMix were smaller than those using the appropriate MCDnCNN model at lower noise levels with smaller orientation angles. For instance, the mean and SD of the AE values derived from Level2 data (orientation angle =0) denoised by the appropriate Delevel2 model were 0.07513 and 0.06069, respectively, whereas those derived from images denoised by DelevelMix were 0.06803 and 0.05413. These values represented a 9.4% lower mean and a 10.8% lower SD of the AE values in using DelevelMix than using Delevel2. However, the situation reversed when the imaging data contained higher levels of noise and/or larger orientation angles, which in reality are more prevail. For example, taking the Level2 data with an orientation angle = π/2, we resulted in the mean and SD of the AE values derived from data denoised by the appropriate Delevel2 model to be 0.06904 and 0.05412, whereas those derived from images denoised by DelevelMix were 0.07350 and 0.05911, respectively. The values based on using Delevel2 were 6.1% (mean) and 8.4% (SD) lower than those based on using DelevelMix, respectively.

When the MWF values calculated based on the Level2 and Level3 data with additionally superimposed noise were compared with the MWF values estimated based on the corresponding noisy Level2 and Level3 data without the added mismatched noise, the maximal mean and SD of the AE of the MWF values were 0.0327 and 0.02555 as long as the mismatched noises were no more than 0.5% (Table 3). However, the results based on an additional 0.75% level of mismatched noise were much worse, although the 0.75% was just slightly higher than 0.5% (Table 3). Similarly, the mean and SD of AE values of the MWF difference between the Level2 data denoised by Delevel2 and the same Level2 data +0.5% denoised by Delevel2 were 0.02639 and 0.02039, whereas those between the Level2 data denoised by Delevel2 and the same Level2 +0.75% data denoised by Delevel2 were much worse at 0.07271 and 0.05656, respectively. A mismatch of ±0.5% noise therefore appeared to be an acceptable limit for using an appropriate denoising model. The results confirmed that the MWF estimates using the appropriate models were good and stable, because all the noise levels employed here approximated to either Level3 or Level2, when the extra noise for testing acceptable levels of mismatched noise was not greater than 0.5%.

In vivo imaging

The data denoised using MCDnCNN models appeared to have helped reduce the disturbance caused by noise in estimating MWF (Figure 4). For example, the statistics based on the raw dataset (noise level =0.0197, corresponding to Level3) showed much larger variations than that based on the corresponding version after denoising (Figure 4). Randomly sampled a voxel from this in vivo image, we see that the T₂* curve of the denoised voxel (the red curve) was very well fitted to the ideal pattern (the blue curve) whereas that based on the raw data without being denoised (green curve) showed much larger variations fluctuating away from the ideal curve (Figure 4A). When the statistics was based on the collection of 30 randomly selected voxels from this in vivo image, we see that the MWF values of their denoised version (the red circles) were lower in both mean and SD than those of the voxels without being denoised (the green circles) (Figure 4B). However, as the noise level was 0.0197, Delevel3 appeared to be the most appropriate denoising model for this particular dataset, and the voxels and the results based on Delevel3 were indeed the smoothest and significantly improved (two sample t-test, P=0.0072), whereas the results generated by Delevel1, Delevel2, Delevel4 and DelevelMix were not as good due to inappropriateness of using the models. In fact, the inappropriate models may lead to either an over- or an under-estimation (Figure 4), as discussed shortly after in the Discussion. The mean and SD of CNR before and after denoising were 45.96±22.56 and 48.76±23.33, respectively, demonstrating no difference (two sample t-test, P=0.7779), which indicated that Delevel3 did not alter the image contrast.

Figure 4 Illustration of the sensitivity of T2* to noise and the effects of the five MCDnCNN models on the T2* decay curve. Panel A shows signal decay as a function of TE at one typical voxel before (green line) and after (red line) denoising. The MWF values (shown in the upper-right rectangular boxes) of the randomly selected voxel were derived from the 5 denoising models (Delevel1, Delevel2, Delevel3, Delevel4 and DelevelMix), respectively. It appears that as the denoising level increased (Delevel1 to Delevel4), the T2* curve became smoother and smoother. Panel B displays the MWF values of 30 randomly selected voxels derived from data before (green circle) and after (red circle) denoising by the 5 models (Delevel1, Delevel2, Delevel3, Delevel4 and DelevelMix), respectively. It appears that as the denoise level increased (Delevel1 to Delevel4), the mean and SD of MWF values became smaller. Note that the overlap of the red and green curves using the DelevelMix model in panel A, and the close matches between the red and green circles in the corresponding figure in panel B (compared with the content in panel B), indicating that the improvement by DelevelMix was indeed minimal. red line, decay curve after denoised; green line, the original decay curve; blue line, complex three-pool model fitted decay curve after denoised; P, P value of two sample t-test between the 30 original and 30 denoised data samples. Delevel, the denoising model corresponding to the noise level; MWF, myelin water fraction; SD, standard deviation; MCDnCNN, multichannel version of the denoising convolutional neural networks.

The comparison of four slices of MWF images based on a randomly picked raw dataset (noise level =0.0197, corresponding to Level3), before and after denoising using the five MCDnCNN models showed that there was substantial random variation in MWF maps derived from the noisy images, as well as artifacts known as “voids” (voxels with very low-MWF values), as indicated by the red rectangle and red arrow in Slice3 and red arrow in Slice2 (Figure 5). After denoising by the MCDnCNN models, the MWF values appeared to have removed most of the artifacts, with those random variations removed and even the voids filled (Figure 5, Slice3 and Slice2, Delevel3). In contrast, the MCDnCNN models were also able to adjust those extraordinarily high MWF values to more appropriate ones (Figure 5, white arrow in Slice3). Moreover, the MWF maps calculated from the MCDnCNN-denoised images exhibited preservation of edges and small structures (Figure 5, green rectangle in Slice3). The results also demonstrated that all the 5 denoising models unanimously reduced random variation although the level of reduction was different. In addition, the MWF values obtained from the noisy images were generally higher than those obtained from the denoised images. The computation times required to denoise the entire dataset using MCDnCNN was about 25 s.

Figure 5 A comparison of four typical slices of MWF maps obtained from an example brain before and after denoising using five MCDnCNN denoising models. The first column to the last column stand for MWF maps at the same 4 slices location of the 6 data groups (the original data set without denoising and the 5 models of Delevel1, Delevel2, Delevel3, Delevel4 and DelevelMix). The two insets beneath Slice3 are the zoomed-in views corresponding to the two rectangular areas outlined in Slice3. Appropriate MCDnCNN models may reduce random variation substantially, fill voids (red rectangle and red arrows), preserve small structural details of the white matter (green rectangle) and correct the overestimated MWF values in some areas (white arrow). The colored bar codes the MWF values from 0% to 35% in the figures. The direction of B0 field was perpendicular to the axial plane (i.e., the view plane in the current case). Delevel, the denoising model corresponding to the noise level; T1W, T1-weughted imaging; MWF, myelin water fraction; MCDnCNN, multichannel version of the denoising convolutional neural networks.

The average noise level of all the 26 subjects were at a mean noise level and SD of 0.01968 and 0.00067 while the average noise level of the 11 subjects for evaluation were at a mean noise level and SD of 0.0202 and 0.000485, which actually corresponded to Level3 (noise level =0.0205). The differences between the MWF values of the five ROIs based on data either denoised by the five respective MCDnCNN models or on data without being denoised were compared (Figure 6). We found that the MWF values derived from data without being denoised were not significantly different from the MWF values obtained from the Delevel1 images except for the forceps major, but were significantly different from the MWF values obtained from the other three models, i.e., Delevel2, Delevel3 and Delevel4 (two-tailed t-test with P<0.05). In addition, the MWF values obtained from the Delevel3 were all lower than those obtained from the noisy images. We compared our MWF results in the ROIs (Table 4) with those in the literatures (10,12,37) and found that for all the ROIs with Delevel3, the results from this current study using the MCDnCNN model were comparable to those in (12) and were the most closely approximated among all the results. Moreover, the SD of the Delevel3 results was the smallest.

Figure 6 In vivo results. Mean and SD MWF values calculated in five different ROIs. The differences of MWF values of the five ROIs between data without denoising and data denoised by the five MCDnCNN models were compared (two-sample t-test). The vertical axis indicates the MWF values, and the horizontal axis indicates the 6 data groups (the original data set without denoising and the 5 models of Delevel1, Delevel2, Delevel3, Delevel4 and DelevelMix). *, P<0.05; **, P<0.001. Delevel, the denoising model corresponding to the noise level; MWF, myelin water fraction; ForcepsMajor, part of Forceps Major; ForcepsMinor, part of Forceps Minor; GCC, part of genu of corpus callosum; IC, part of internal capsules; SCC, part of splenium of corpus callosum; SD, standard deviation; ROIs, regions-of-interests; MCDnCNN, multichannel version of the denoising convolutional neural networks.

Discussion

Estimating MWF is a challenging work, as it involves dealing with both noisy imaging data and complex data models. MCDnCNN appears to be able to address the challenge quite well, which may significantly enhance the quality of deriving MWF maps from mGRE imaging data. Our numerical analysis demonstrated markedly reduced AE in MWF maps calculated based on data denoised by an appropriate MCDnCNN model as compared to noised images (Figure 3, Tables 1,2). Furthermore, analysis using in vivo data showed that estimation of MWF was markedly improved through use of the appropriate MCDnCNN models (Figures 4-6). Indeed, our results indicated that MWF using MCDnCNN benefits the estimation of MWF maps in preserving anatomic details and in delineation of local myelin patterns as compared to results obtained directly from noisy imaging data.

The general applicability of the MCDnCNN models has been validated using different MRI datasets (19), that MCDnCNN is able to estimate any additional noise (e.g., Rician noise) superimposed onto the original image. In the subsequent process of denoising, the original image will be treated as a clean image, even if it appears noisy. Previous papers reported that a MCDnCNN could denoise different data bearing the same noise distribution (19,22). The general applicability is correlated to residual learning formulation and BN, which are two key features of MCDnCNN differentiating it from other types of deep leaning networks (19,22). The residual learning formulation is adopted to learn a residual mapping with the operations in the hidden layers but not to learn a mapping function to predict the latent clean image (20,38,39). This makes the inputs to each next layer specific-noise-like distributed, so that the inputs to the next layers are less and less correlated with the latent image content (20). For instance, in Rician denoising the residual image is associated with the Rician distribution instead of the distribution of the latent clean images. Although the contrast and noise level change significantly depend on TE, the noise distribution is typically stable across the entire decay along all the TEs (40). After training, the network would be able to estimate the noise in each TE image. Thus, the contrast and noise level would have a less impact on denoising while the denoising model is appropriate. Therefore, MCDnCNN is suitable for denoising the images at each TE of the mGRE datasets. Using residual learning, we aimed to learn the noise distribution for obtaining the latent clean images but not directly the distribution of the latent clean images. However, we should try other powerful networks in the future, such as the U-net (41), to predict the distribution of latent clean images directly to make a comparison with our current work. In addition, training the model using clean mGRE data is certainly preferred, whenever the scanning technology of acquiring clean mGRE data becomes available. We expect a comparison based on such a model can be done with our work in the near future.

Simulations were conducted using a hollow cylinder model that incorporated susceptibility anisotropy fiber orientation dependencies. We have seen that appropriate MCDnCNN models generated smaller AE in MWF maps than inappropriate MCDnCNN models did (Figure 3, Tables 1,2). Moreover, MWF estimation is sensitive to noise and denoising levels. MWF values derived from data treated by inappropriate MCDnCNN models were always overestimated or underestimated. MWF values derived from data without being denoised had the same issue (Figure 6). The MWF values at voxels within the five ROIs demonstrated significant difference between the data with and without denoising. This may be attributed to the complexity of fitting the 3-pool model with ten unknown variables to a decay curve based on the mGRE data of 24 echoes, however, within only a limited range of TEs, which is actually an inherently ill-conditioned problem (42). Consequently, the calculation may frequently fall into local minima and thus became sensitive to noise in fitting the curve. However, the results in the experiment based on simulation data denoised by the appropriate MCDnCNN models met our expectation that the MWF values should fall within 5–18%. Moreover, the results in the experiments based on in vivo data denoised by the appropriate MCDnCNN models were generally in good agreement with those previously reported in the literature (Table 4). The observation showed that using the appropriate denoising models may effectively alleviate the ill-conditioned problem. In fact, MCDnCNN was designed to predict the residual images (noise), i.e., the difference between the noisy observation and the latent clean image (20), thus adopting MCDnCNN was expected to remove the noise from the input data to avoid the aforementioned ill-conditioned problem. In other words, MCDnCNN may learn the pattern of the noise from the input data. The noise of the same pattern can therefore be removed from the mGRE data at each TE map, yet, without decreasing the resolution of the images.

The MWF values derived from the denoised in vivo data treated by the MCDnCNN models agreed well with the results derived from the data in the simulation experiments, whenever an appropriate denoising model was applied in accordance with the noise level of the corresponding dataset. It infers that to avoid overestimating or underestimating MWF in the in vivo data, we must first estimate the noise level of the data, and then particularly choose the Delevel model that is appropriate for the noise level of the data. In Figure 5, we see that using an appropriate MCDnCNN model (i.e., Delevel3) may improve the quality of MWF mapping, including filling the voids, reducing over-estimation and preserving structural details. Although the void (red rectangle in Slice3; Figure 5) was at the boundary of white matter adjacent to neighboring gray matter, which is visible better in T1-weighted image of the same subject, myelin component was expected. Similar but more convincing evidence is presented at the tips of the red arrows in Figure 5 on both Slice2 and Slice3. The areas pointed by the white arrow in Slice3 contained lots of myelinated fiber bundles. However, most of these fibers were approximately parallel or perpendicular to the B0 direction in our case (B0 was perpendicular to the axial plane in Figure 5). Previous studies reported that the angle between the myelination fiber bundles and the B0 direction will affect the accuracy of estimating MWF (9,10,12), just like what we have seen in our experiments using the simulated data. In line with the findings in these papers, the myelin components in the fiber bundles that were either parallel or perpendicular to the B0 direction would achieve lower SD in estimating the MWF values. According to the latest research, the MWF value range in this portion of the fiber bundles in the brain is about 9–15% (43), and our results of the data before denoising were obviously higher while the results of the data after denoising by Delevel3 were consistent with the reported in the literature. Therefore, we reasonably believe that our results showed an improvement of accuracy in MWF estimation.

In contrast, the other four denoising models (Delevel1, Delevel2, Delevel4 and DelevelMix) performed poorly. For example, the same data denoised by the four models failed to estimate meaningful MWF values to fill the voids (Figure 5, areas at the red arrows in both Slice2 and Slice3). And in some areas, data denoised by Delevel1, Delevel2 and DelevelMix resulted in MWF values that were obviously too higher than the reported range (Figure 5, areas indicated by the white arrows), while data denoised by Delevel4 resulted in MWF values that were obviously lower than the reported range (Figure 5). Moreover, extreme results appeared in the ROI analysis of the data denoised by Delevel1 and DelevelMix that, for example, the mean of MWF values in the forceps minor were overly higher than those derived from original data, which were already higher than the reported results (Table 4) (12,35). Thus, the inappropriate models led to denoised data that turned out either an over- or an under-estimation of MWF values.

Checking against the corresponding T1-weighted image (last row, Figure 5) in which the structural boundaries are more intuitive, we may see that the small structural details in white matter were well preserved in the data denoised by Delevel3 (Figure 5, green rectangle in Slice3). Yet, we noticed that some structures, such as the genu in the frontal lobe, disappeared in the results based on the data denoised by the MCDnCNN models. The disappearance here was also observed in data of high SNR acquired using an ultra-high field strength at 7 T (10). In fact, investigators speculated that such disappearance is not caused by denoising but is related to B0 inhomogeneity (8,17). Our main purpose in the current work was to study the benefits of employing appropriately denoised data in MWF estimation, however, denoising cannot correct artifacts caused by B0 inhomogeneity. The experiments have shown that the mean values of the calculated MWF values based on data denoised using the appropriate models were generally in good agreement with those previously reported in the literature (Table 4). However, the SD values of the MWF, especially in the genu of corpus callosum, were slightly higher than those reported, although they were still within the acceptable range according to Alonso-Ortiz’s report (44). We speculate that this phenomenon was due to a higher spatial resolution, and thereby a reduced SNR and inherently the aforementioned issue of ill-conditioned problem. Further experiments and analyses would be necessary to double check this speculation.

We have seen that the MCDnCNN models may permit substantially improved quality of MWF maps derived from mGRE imaging datasets. Employing MCDnCNN models has a number of advantages. First of all, significantly improved computational efficiency: MCDnCNN took only about 25 s to denoise an entire series of brain volumes, whereas few minutes were typically needed in denoising such a series of brain volumes using conventional methods (19,45,46). Second, quality improvement: employing the appropriate MDCnCNN models may well preserve brain details and image resolution. Third, convenience in training this network (19): the implementation of the network and its training are technically straightforward, which rely on only one single user-defined parameter, i.e., the noise level, defined as the noise at the last TE divided by the mean signal at the first TE. However, an appropriate MDCnCNN model suitable for this noise level must be used in this scenario, because using an inappropriate model will lead to an overestimation or underestimation, as demonstrated in the experiments. Previous papers smoothed the images using a non-local means filter (47) or an anisotropic diffusion filter (ADF), but the denoising parameter was determined by the SD of the noise in the 3D volumetric image acquired at the first echo (48). However, in the denoising practice, we noted that the T₂* decay curve typically had larger deviation in the last TE images (Figure 4), making the data of the last TE image relatively more seriously biased. To effectively control and reduce such deviation and bias, and subsequently to maximize the tolerance of noise deviation in different TE images, we used only first TE and the last TE images to estimate the noise level. Our experiments indicated that whenever the denoised level (using an appropriate denoising model) corresponded to the noise level, credible results may be achieved. Although estimating the noise level in advance imposes an additional step, the slight overhead will not offset the huge advantages of employing MCDnCNN. Moreover, the SNR increase without sacrificing image resolution in the denoised images is desired in clinical applications, known as enhanced clinical applicability of advanced signal models exhibiting high sensitivity to noise (49,50).

It should be noted that the generalized MCDnCNN models, i.e., DelevelMix, showed better performance than an appropriate MDCnCNN model at lower noise level(s) with lower orientation angles in the experiments using the simulated data (Tables 1,2). However, when in vivo data were used, DelevelMix showed no advantages at all over Delevel3 (Table 4 and Figures 4-6). These results were in line with those previously reported (19,20). Although the generalized convolutional neural networks may satisfy specific denoising tasks with certain levels of noise, their performance is generally not as stable as that of the noise-specific networks under the same conditions. Experiment results based on simulated and in vivo data showed that 3CCT₂* is sensitive to noise that even a small disturbance may bias the results (Figure 4 and Table 3). However, our experiments also have provided helpful hints for training of the generalized MCDnCNN models, such as the sensitivity of the angle between the running orientation of the underlying fibers and B0 direction, which should be taken care of more carefully in the analyses.

The sequence of mGRE is being widely used in estimating MWF because of the advantages of the T₂* MWI model, including the availability of the intercorrelation between echoes for use (7,8). However, the challenge turns to be the issue of SNR, due to the reduced size of the voxels when mGRE is applied for acquisition of images with a high spatial resolution. The low SNR will lead to the MWF estimation more prone to systematic errors using the direct complex fitting procedure, and consequently severely compromising stability and reliability of the estimation (12). Data with intrinsically low SNR based on small-sized voxels used in such correction procedures (17,51) may make the estimation highly unreliable. Although the method of using QSM and Eq. [2] in the current work to estimate tissue frequency and phase is a simplified approach, the signal decay accounted for QSM-derived phases may partly correct errors induced by magnetic field inhomogeneity to certain extent and thus may improve the estimation of the local phase introduced by the white matter tissue, thereby making 3CCT₂* a better model in deriving MWF maps (12,33,52). As the focus of the current work was to examine the ability of MCDnCNN in estimating MWF, we did not further scrutinize this aspect but considered leaving it to a separated work in the near future. Although our procedure did not completely resolve some inherent QSM problems, such as motion artifacts caused by unbalanced gradient and pulsation (31,53,54), a recent paper succeeded in training a deep convolutional neural network to reduce the motion artifact in QSM (55). In addition, our proposed method only affected the magnitude images but not the complex images, therefore had little influence on phase images. A deep convolutional neural network trained to denoise the complex images was proposed in a recent paper (56), which might improve the phase images. Such networks aiming to overcome the remaining QSM problems and to improve the quality of both the magnitude images and the phase images shall be tested in the future work. Also, previous studies had used mGRE to obtain susceptibility maps by employing a pipeline involving QSM (12,33,52,57). Wu et al. and Sood et al. (12,33) demonstrated that the accuracy in MWF estimation can get improved using the QSM pipeline with the STI suite by overcoming non-local phase contrast and orientation dependence, along with denoising the phase component. We therefore followed the same schema as an initial step of the effort, though it was relatively straightforward.

Some limitations existed in our study. First, all the involved subjects were healthy adults and therefore we did not have a chance to test its clinical usefulness using patient data, although the results have demonstrated the effectiveness of the proposed method. Future studies shall try validating the proposed method on patient data such as those of multiple sclerosis, schizophrenia, Alzheimer disease and stroke. Second, the simplified 1D representation of the hollow cylinder model (21) could introduce certain error into the estimation of T₂*, especially for MW components in fibers that run in an orientation with a relatively large angle relative to the direction of the main magnetic field B0. Such error, if any, could eventually propagate to the MWF maps (43). For instance, SD values of the AE between the calculated MWF and the ground truth MWF close to the mean values could imply a wide distribution of the calculated MWF. These phenomena were in line with previous results that 3CCT₂* had a lager variability (10). Although we have demonstrated through experiments with simulation data that the appropriate MCDnCNN denoising model may achieve ideal results, we shall further improve the model to cater better for the realistic conditions, and thereby address the possible inaccurate estimation of T₂* and the consequent bias in the MWF maps (58,59). Nevertheless, past studies have pointed out that it would mostly affect a magnitude fitting model at higher field strengths (10,43). Third, the density of noise levels used for training may need to be increased. The models we trained were at Rician noise levels with a 2% increment (19). This incremental step seems to be a little too large for high-resolution images though. In our future study, we shall test training the denoising models using a smaller increment, such as using 1%, 0.5%, and 0.1%, which would allow us to study the denoising effects with a finer observation. Fourth, we used Rician noise to simplify a non-central chi distribution instead of using other more complex noise models while using GRAPPA to reconstruct the magnitude images, considering that GRAPPA reconstruction was non-stationary (60), the noise distribution on mGRE was complex, and our major purpose was to improve MWF calculation using MCDnCNN. Our future efforts should consider using a more complicated noise model. Another limitation was the limited range of TE permitted by our scanner, which affected the visualization and also possibly the degree of accurate estimation of the multi-exponential shape of the decay curve (Figure 4). While a previous study (43) showed that the MWF value stabilizes at TE around 40 ms using the complex three-pool model, we believe that an extended TE range would benefit the study as information will be more comprehensive. We should therefore expand the TE value range, once the scanning technology advances may allow us to do so.

Conclusions

In conclusion, the accuracy and quality of in vivo estimation of MWF in the human brain from high spatial resolution with limited TEs were markedly improved through use of the appropriate MCDnCNN denoising models. The use of MCDnCNN denoising models may contribute significantly to the goal of obtaining high-quality MWF mapping efficiently that is useful in clinical practice and settings.

Acknowledgments

Funding: This work was supported in part by a grant by China National Key R&D Program (2017YFC1308500/2017YFC1308502) and also grants sponsored by China National Science Foundation (81471651, 11835003, 81471734).

Footnote

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://dx.doi.org/10.21037/qims-21-404). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). This study was approved by the local Ethics Review Board of East China Normal University (No. HR227-2019) and informed consent was obtained from all the participants.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Bakshi R, Thompson AJ, Rocca MA, Pelletier D, Dousset V, Barkhof F, Inglese M, Guttmann CR, Horsfield MA, Filippi M. MRI in multiple sclerosis: current status and future prospects. Lancet Neurol 2008;7:615-25. [Crossref] [PubMed]
2. Flynn SW, Lang DJ, Mackay AL, Goghari V, Vavasour IM, Whittall KP, Smith GN, Arango V, Mann JJ, Dwork AJ, Falkai P, Honer WG. Abnormalities of myelination in schizophrenia detected in vivo with MRI, and post-mortem with analysis of oligodendrocyte proteins. Mol Psychiatry 2003;8:811-20. [Crossref] [PubMed]
3. Bouhrara M, Reiter DA, Bergeron CM, Zukley LM, Ferrucci L, Resnick SM, Spencer RG. Evidence of demyelination in mild cognitive impairment and dementia using a direct and specific magnetic resonance imaging measure of myelin content. Alzheimers Dement 2018;14:998-1004. [Crossref] [PubMed]
4. Borich MR, Mackay AL, Vavasour IM, Rauscher A, Boyd LA. Evaluation of white matter myelin water fraction in chronic stroke. Neuroimage Clin 2013;2:569-80. [Crossref] [PubMed]
5. Ma YJ, Jang H, Chang EY, Hiniker A, Head BP, Lee RR, Corey-Bloom J, Bydder GM, Du J. Ultrashort echo time (UTE) magnetic resonance imaging of myelin: technical developments and challenges. Quant Imaging Med Surg 2020;10:1186-203. [Crossref] [PubMed]
6. MacKay A, Whittall K, Adler J, Li D, Paty D, Graeb D. In vivo visualization of myelin water in brain by magnetic resonance. Magn Reson Med 1994;31:673-7. [Crossref] [PubMed]
7. Du YP, Chu R, Hwang D, Brown MS, Kleinschmidt-DeMasters BK, Singel D, Simon JH. 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:865-70. [Crossref] [PubMed]
8. Hwang D, Kim DH, Du YP. In vivo multi-slice mapping of myelin water content using T2* decay. Neuroimage 2010;52:198-204. [Crossref] [PubMed]
9. Nam Y, Lee J, Hwang D, Kim DH. Improved estimation of myelin water fraction using complex model fitting. Neuroimage 2015;116:214-21. [Crossref] [PubMed]
10. Alonso-Ortiz E, Levesque IR, Pike GB. Impact of magnetic susceptibility anisotropy at 3 T and 7 T on T2*-based myelin water fraction imaging. Neuroimage 2018;182:370-8. [Crossref] [PubMed]
11. Lee J, Nam Y, Choi JY, Kim EY, Oh SH, Kim DH. Mechanisms of T2 * anisotropy and gradient echo myelin water imaging. NMR Biomed 2017;30: [Crossref] [PubMed]
12. Wu Z, He H, Sun Y, Du Y, Zhong J. High resolution myelin water imaging incorporating local tissue susceptibility analysis. Magn Reson Imaging 2017;42:107-13. [Crossref] [PubMed]
13. Istratov AA, Vyvenko OF. Exponential analysis in physical phenomena. Rev Sci Instrum 1999;70:1233-57. [Crossref]
14. Alonso-Ortiz E, Levesque IR, Pike GB. MRI-based myelin water imaging: A technical review. Magn Reson Med 2015;73:70-81. [Crossref] [PubMed]
15. Hwang D, Chung H, Nam Y, Du YP, Jang U. Robust mapping of the myelin water fraction in the presence of noise: synergic combination of anisotropic diffusion filter and spatially regularized nonnegative least squares algorithm. J Magn Reson Imaging 2011;34:189-95. [Crossref] [PubMed]
16. Jones CK, Whittall KP, MacKay AL. Robust myelin water quantification: averaging vs. spatial filtering. Magn Reson Med 2003;50:206-9. [Crossref] [PubMed]
17. Alonso-Ortiz E, Levesque IR, Paquin R, Pike GB. Field inhomogeneity correction for gradient echo myelin water fraction imaging. Magn Reson Med 2017;78:49-57. [Crossref] [PubMed]
18. Bouhrara M, Bonny JM, Ashinsky BG, Maring MC, Spencer RG. Noise Estimation and Reduction in Magnetic Resonance Imaging Using a New Multispectral Nonlocal Maximum-likelihood Filter. IEEE Trans Med Imaging 2017;36:181-93. [Crossref] [PubMed]
19. Jiang D, Dou W, Vosters L, Xu X, Sun Y, Tan T. Denoising of 3D magnetic resonance images with multi-channel residual learning of convolutional neural network. Jpn J Radiol 2018;36:566-74. [Crossref] [PubMed]
20. Zhang K, Zuo W, Chen Y, Meng D, Zhang L. Beyond a Gaussian Denoiser: Residual Learning of Deep CNN for Image Denoising. IEEE Trans Image Process 2017;26:3142-55. [Crossref] [PubMed]
21. Wharton S, Bowtell R. Fiber orientation-dependent white matter contrast in gradient echo MRI. Proc Natl Acad Sci U S A 2012;109:18559-64. [Crossref] [PubMed]
22. Li W, Wang N, Yu F, Han H, Cao W, Romero R, Tantiwongkosi B, Duong TQ, Liu C. A method for estimating and removing streaking artifacts in quantitative susceptibility mapping. Neuroimage 2015;108:111-22. [Crossref] [PubMed]
23. He K, Zhang X, Ren S, Sun J. Deep Residual Learning for Image Recognition. IEEE Conference on Computer Vision and Pattern Recognition 2016:770-8.
24. Krizhevsky A, Sutskever I, Hinton GE. Imagenet classification with deep convolutional neural networks. Adv Neural Inf Process Syst 2012;25:1097-105.
25. Ioffe S, Szegedy C. Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift. JMLRorg 2015.
26. Boureau YL, Ponce J, Lecun Y. editors. A Theoretical Analysis of Feature Pooling in Visual Recognition. ICML Proceedings, 27th International Conference on Machine Learning, Haifa: ICML, 2010:111-8.
27. Kingma DP, Ba J. Adam: A method for stochastic optimization. arXiv 2014;arXiv:14126980.
28. Cotter A, Shamir O, Srebro N, Sridharan K. Better mini-batch algorithms via accelerated gradient methods. arXiv 2011;arXiv:11064574.
29. Vedaldi A, Lux M, Bertini M. MatConvNet: CNNs are also for MATLAB users. ACM Sigmultimedia Records 2018;10:9. [Crossref]
30. Reichenbach JR, Venkatesan R, Yablonskiy DA, Thompson MR, Lai S, Haacke EM. Theory and application of static field inhomogeneity effects in gradient-echo imaging. J Magn Reson Imaging 1997;7:266-79. [Crossref] [PubMed]
31. Schweser F, Deistung A, Lehr BW, Reichenbach JR. Quantitative imaging of intrinsic magnetic tissue properties using MRI signal phase: an approach to in vivo brain iron metabolism? Neuroimage 2011;54:2789-807. [Crossref] [PubMed]
32. Schweser F, Deistung A, Reichenbach JR. Foundations of MRI phase imaging and processing for Quantitative Susceptibility Mapping (QSM). Z Med Phys 2016;26:6-34. [Crossref] [PubMed]
33. Sood S, Urriola J, Reutens D, O'Brien K, Bollmann S, Barth M, Vegh V. Echo time-dependent quantitative susceptibility mapping contains information on tissue properties. Magn Reson Med 2017;77:1946-58. [Crossref] [PubMed]
34. Sati P, van Gelderen P, Silva AC, Reich DS, Merkle H, de Zwart JA, Duyn JH. Micro-compartment specific T2* relaxation in the brain. Neuroimage 2013;77:268-78. [Crossref] [PubMed]
35. Kim D, Lee HM, Oh SH, Lee J. Probing signal phase in direct visualization of short transverse relaxation time component (ViSTa). Magn Reson Med 2015;74:499-505. [Crossref] [PubMed]
36. Magnotta VA, Friedman L. FIRST BIRN. Measurement of Signal-to-Noise and Contrast-to-Noise in the fBIRN Multicenter Imaging Study. J Digit Imaging 2006;19:140-7. [Crossref] [PubMed]
37. Whittall KP, MacKay AL, Graeb DA, Nugent RA, Li DK, Paty DW. In vivo measurement of T2 distributions and water contents in normal human brain. Magn Reson Med 1997;37:34-43. [Crossref] [PubMed]
38. Schmidt U, Roth S. Shrinkage Fields for Effective Image Restoration. IEEE Conference on Computer Vision and Pattern Recognition 2014:2774-81.
39. Burger HC, Schuler CJ, Harmeling S. editors. Image denoising: Can plain neural networks compete with BM3D? Computer Vision and Pattern Recognition (CVPR), Providence, RI, USA: IEEE, 2012: 2392-9.
40. Nam Y, Kim DH, Lee J. Physiological noise compensation in gradient-echo myelin water imaging. Neuroimage 2015;120:345-9. [Crossref] [PubMed]
41. Batson J, Royer L. Noise2self: Blind denoising by self-supervision. PMLR 2019;97:524-33.
42. van Gelderen P, de Zwart JA, Lee J, Sati P, Reich DS, Duyn JH. Nonexponential T₂ decay in white matter. Magn Reson Med 2012;67:110-7. [Crossref] [PubMed]
43. Lee H, Nam Y, Kim DH. Echo time-range effects on gradient-echo based myelin water fraction mapping at 3T. Magn Reson Med 2019;81:2799-807. [Crossref] [PubMed]
44. Alonso-Ortiz E, Levesque IR, Pike GB. Multi-gradient-echo myelin water fraction imaging: Comparison to the multi-echo-spin-echo technique. Magn Reson Med 2018;79:1439-46. [Crossref] [PubMed]
45. Manjón JV, Coupé P, Buades A, Louis Collins D, Robles M. New methods for MRI denoising based on sparseness and self-similarity. Med Image Anal 2012;16:18-27. [Crossref] [PubMed]
46. Wu X, Yang Z, Peng J, Zhou J. Global denoising for 3D MRI. Biomed Eng Online 2016;15:54. [Crossref] [PubMed]
47. Coupe P, Yger P, Prima S, Hellier P, Kervrann C, Barillot C. An optimized blockwise nonlocal means denoising filter for 3-D magnetic resonance images. IEEE Trans Med Imaging 2008;27:425-41. [Crossref] [PubMed]
48. Gerig G, Kubler O, Kikinis R, Jolesz FA. Nonlinear anisotropic filtering of MRI data. IEEE Trans Med Imaging 1992;11:221-32. [Crossref] [PubMed]
49. Avram AV, Sarlls JE, Barnett AS, Özarslan E, Thomas C, Irfanoglu MO, Hutchinson E, Pierpaoli C, Basser PJ. Clinical feasibility of using mean apparent propagator (MAP) MRI to characterize brain tissue microstructure. Neuroimage 2016;127:422-34. [Crossref] [PubMed]
50. Steven AJ, Zhuo J, Melhem ER. Diffusion kurtosis imaging: an emerging technique for evaluating the microstructural environment of the brain. AJR Am J Roentgenol 2014;202:W26-33. [Crossref] [PubMed]
51. Lee H, Nam Y, Lee HJ, Hsu JJ, Henry RG, Kim DH. Improved three-dimensional multi-echo gradient echo based myelin water fraction mapping with phase related artifact correction. Neuroimage 2018;169:1-10. [Crossref] [PubMed]
52. Kadamangudi S, Reutens D, Sood S, Vegh V. Signal compartments in ultra-high field multi-echo gradient echo MRI reflect underlying tissue microstructure in the brain. Neuroimage 2018;178:403-13. [Crossref] [PubMed]
53. Fortier V, Levesque IR. Phase processing for quantitative susceptibility mapping of regions with large susceptibility and lack of signal. Magn Reson Med 2018;79:3103-13. [Crossref] [PubMed]
54. Li W, Wu B, Liu C. Quantitative susceptibility mapping of human brain reflects spatial variation in tissue composition. Neuroimage 2011;55:1645-56. [Crossref] [PubMed]
55. Li C, Zhang H, Zhang J, Spincemaille P, Nguyen TD, Wang Y. Motion Artifact Reduction in Quantitative Susceptibility Mapping using Deep Neural Network. arXiv 2021;arXiv:210501746.
56. Quan Y, Chen Y, Shao Y, Teng H, Xu Y, Ji H. Image denoising using complex-valued deep CNN. Pattern Recognition 2021;111:107639. [Crossref]
57. Nissi MJ, Tóth F, Wang L, Carlson CS, Ellermann JM. Improved Visualization of Cartilage Canals Using Quantitative Susceptibility Mapping. PLoS One 2015;10:e0132167. [Crossref] [PubMed]
58. Xu T, Foxley S, Kleinnijenhuis M, Chen WC, Miller KL. The effect of realistic geometries on the susceptibility-weighted MR signal in white matter. Magn Reson Med 2018;79:489-500. [Crossref] [PubMed]
59. Yablonskiy DA, Sukstanskii AL. Biophysical mechanisms of myelin-induced water frequency shifts. Magn Reson Med 2014;71:1956-8. [Crossref] [PubMed]
60. Aja-Fernández S, Tristán-Vega A, Hoge WS. Statistical noise analysis in GRAPPA using a parametrized noncentral Chi approximation model. Magn Reson Med 2011;65:1195-206. [Crossref] [PubMed]