The Utility of a Convolutional Neural Network for Generating a Myelin Volume Index Map from Rapid Simultaneous Relaxometry Imaging.

PURPOSE: A current algorithm to obtain a synthetic myelin volume fraction map (SyMVF) from rapid simultaneous relaxometry imaging (RSRI) has a potential problem, that it does not incorporate information from surrounding pixels. The purpose of this study was to develop a method that utilizes a convolutional neural network (CNN) to overcome this problem.
METHODS: RSRI and magnetization transfer images from 20 healthy volunteers were included. A CNN was trained to reconstruct RSRI-related metric maps into a myelin volume-related index (generated myelin volume index: GenMVI) map using the MVI map calculated from magnetization transfer images (MTMVI) as reference. The SyMVF and GenMVI maps were statistically compared by testing how well they correlated with the MTMVI map. The correlations were evaluated based on: (i) averaged values obtained from 164 atlas-based ROIs, and (ii) pixel-based comparison for ROIs defined in four different tissue types (cortical and subcortical gray matter, white matter, and whole brain).
RESULTS: For atlas-based ROIs, the overall correlation with the MTMVI map was higher for the GenMVI map than for the SyMVF map. In the pixel-based comparison, correlation with the MTMVI map was stronger for the GenMVI map than for the SyMVF map, and the difference in the distribution for the volunteers was significant (Wilcoxon sign-rank test, P < 0.001) in all tissue types.
CONCLUSION: The proposed method is useful, as it can incorporate more specific information about local tissue properties than the existing method. However, clinical validation is necessary.

Introduction

Myelin is an important component of central nervous system. It maintains the integrity of the neuron fibers and accelerates the propagation of action potentials, which are inevitable to maintain the appropriate function.[1,2] Thus, measuring myelin is important to diagnose or monitor the disorders and diseases that results from myelin dysfunction. MRI is the most important tool to measure myelin in vivo. Its usefulness is already established in evaluating the development and aging of humans and for assessing the progression of degenerative or demyelinating diseases.[3-5] However, the scan time is one of the important issues that might become problematic for general clinical practice.[3] Conventional imaging sequences (e.g. T1- and T2-weighted image) are very sensitive to tissue contrast, but not as well for the detailed tissue properties such as myelin content. So, additional scan is required in some of the major techniques to obtain myelin volume index (MVI) including the one based on magnetization-saturation imaging (MTsat).[6,7]

Recently, rapid simultaneous relaxometry imaging (RSRI) was developed to measure the longitudinal relaxation rate (R1), transverse relaxation rate (R2), proton density (PD), and the local B1 field from a single scan within an acceptable scan-time.[8] Here, the metrics can be used to estimate the myelin volume fraction (MVF) using a commercial software package SyMRI (Synthetic MR, Linköping, Sweden).[9] This function has a great potential merit for myelin assessment in the general clinic because it does not require additional scan time to provide the myelin volume estimate. In addition, the MVF estimated with SyMRI [i.e. synthetic myelin volume fraction (SyMVF)] was reported as being highly correlated with MTsat-based myelin measurement, namely the magnetization transfer-based MVI (MTMVI).[3] Furthermore, the usefulness of SyMVF has already been noted in several applications.[10-13] However, there might be room for improvement in the algorithm used to estimate a SyMVF map. The SyMVF metric is determined pixel-by-pixel from a lookup-table that connects combinations of R1, R2, and PD values to a MVF,[11] meaning that a pixel with a certain combination of R1, R2, and PD values is always assigned the same MVF without considering any local properties. As tissue structure differs in different areas of the brain (e.g. neuron count,[14] neuron fiber radius,[14,15] iron deposition), this could lead to inaccuracy when generating the SyMVF map. Adding information about local tissue properties may strengthen the accuracy and reliability of the output MVF map.

Recently, the convolutional neural network (CNN) technique achieved great success for image segmentation of many areas in the human body.[16-19] As the shape of the data of each layer of a CNN is generally unrestricted, taking the processing stream from one CNN and adding it to the function of another CNN is possible (e.g. Litjens et al.,[17] Çiçek et al.,[18] and Perone et al.[19]). Based on this idea, this study combines a CNN for segmentation with another simple CNN designed for general non-linear reconstruction. The purpose to use the first CNN for segmentation is to extract the information of the local structure from the R1, R2, and PD maps. The second simple CNN was used to combine the local structure information that was provided from the first CNN with the existing SyMVF map to generate a new MVI (GenMVI) map. From this architecture, the GenMVI map is expected that it can incorporate the information of the local tissue and become more specific to the characteristics of the tissue in each pixel than SyMVF.

The purpose of this study was to evaluate the usefulness of this method for estimating myelin-volume in human brain.

Materials and Methods

This study was performed as a part of a prospective study. The data used in this study was also used in another previous study.[3] This study was approved by the IRB of Juntendo University. Written informed consent was obtained from all participants.

Study participants

Twenty healthy volunteers those went through the scans explained in the next section for a previous study[3] were included as candidates for the subject of this study. The volunteers were, nine males (25–67 years, mean 53.2 years) and 11 females (44–71 years, mean 57.0 years), without neurological or psychological history. The images acquired from the volunteers were screened by two board-certified radiologists (Y.T. and A.H., 12 and 5 years of experience interpreting brain MRI, respectively) to exclude the volunteers with moderate-to-severe white-matter ischemic lesions (Fazekas grade 2 or more[20]), asymptomatic cerebral infarction, or regional brain atrophy. No volunteer was excluded from the subject.

Image acquisition and data processing to generate SyMVF and MTMVI maps

All scans in this study were performed by 3T MRI scanner (MAGNETOM Prisma, Siemens Healthcare, Erlangen, Ger-many) using a 64-channnel head coil. Images were acquired using the imaging sequence of QRAPMASTER (quantification of relaxation times and proton density by multi-echo acquisition of saturation-recovery using turbo spin-echo readout),[3] which has two different TEs (i.e. 22 and 99 ms) and four different saturation delay times (i.e. 170, 620, 1970, and 4220 ms) in a single scan. The other major parameters for QRAPMASTER were: TR 4250 ms; FOV 230 × 186 mm; matrix 320 × 186; slice thickness/gap 4.0/1.0 mm. The acquired images were processed using SyMRI 8.0 software (SyntheticMR, Linkoping, Sweden) to obtain R1, R2, PD, and SyMVF maps. A brain-area probability (BAP) map and synthetic T1-weighted image were also automatically generated during this process.

Imaging with magnetization-transfer (MT) weighting was performed as a first step to obtain MTMVI images. Images were acquired using FLASH sequence with T1-, PD-, and magnetization transfer-weightings. TR and excitation flip angle were set at 10 ms and 13° for T1-weighted images, and 24 ms and 4° for PD- and MT-weighted images. An off-resonance Gaussian-shaped RF pulse (frequency offset from water resonance 1.2 kHz, pulse duration 9.984 ms, and nominal flip angle 50°) was adopted for the MT-weighted images. The other major parameters were: FOV 224 × 224 mm; matrix 128 × 128; slice thickness 1.8 mm. The MTsat map was calculated from the images as described in a previous report,[3] and it was then scaled as described in the next two sections to generate a final MTMVI map.

ROI Definition

The Johns Hopkins University ICBM-DTI-81 WM labels atlas[21,22] and Automated Anatomical Labeling atlas[23,24] were used to define 48 local ROIs for the WM area, and 108 and eight ROIs for the cortical and subcortical GM areas, respectively.

Local ROIs of the atlases were registered to SyMVF volumes for each volunteer. First, the synthetic T1-weighted image volume of a volunteer was registered to the MNI152 template using the FMRIB Software Library (FSL) linear image registration tools and non-linear image registration tools.[25,26] The warp function defined in this registration was then inverted to warp the atlas ROIs (total 164 local ROIs) to fit the volunteer’s SyMVF space. In addition, the warped ROIs were grouped and merged to form another set of ROIs: cortical GM (ROIcGM), subcortical GM (ROIsGM) and WM (ROIWM). These three ROIs were eroded once with an eight-connected-neighborhood rule to avoid partial-volume effects at the margins of each tissue type. Note that when using one of the 164 local ROIs in the analysis, pixels within the ROI that had been eroded from either ROIcGM, ROIsGM, or ROIWM were removed for the same reason. In addition, a whole brain ROI (ROIWB) was created by merging ROIcGM, ROIsGM, and ROIWM. Furthermore, local ROIs corresponding to the genu, body, and splenium of the corpus-callosum were merged to form a single ROI for the corpus-callosum (ROICC).

ROI definition for the MTsat volume was performed in the same way as for the SyMVF map, except that a 3D T1-weighted image volume was used instead of a synthetic T1-weighted image volume.

MTMVI maps

Magnetization-saturation imaging volumes were scaled to create the MTMVI map so that all myelin-related map images in this study are on the same scale. For each MTsat and SyMVF volume pair, pixels in the WM area were extracted and averaged using the ROIWM defined for each volunteer in the previous section. The MTMVI volume was created by multiplying MTsat volume by a constant scaling value so that the average values of pixels in ROIWM were equal for the SyMVF and MTMVI volumes.

Finally, each MTMVI volume was nonlinearly registered to the SyMVF volume space of the same volunteer. The Advanced Normalization Tools (antsRegistrationSyNQuick.sh, http://stnava.github.io/ANTs/) package[27] was used for this purpose.

Deep-learning-based method to obtain myelin volume index

Image preparation

Computing procedures described in this section were performed using our in-house software running on MATLAB 2017b® (Mathworks, Natick, USA). To create a dataset for training, first, 32 × 32-pixel patch images were randomly subsampled from each slice of all image volumes (i.e. R1, R2, PD, BAP, SyMVF, and MTMVI) corresponding to each volunteer. Patches including brain area (defined from BAP map as pixels with P > 0.95) of less than half of the whole area were excluded. Finally, approximately 6000 patch sets were subsampled from each volunteer. All subsampled patches were resized to 128 × 128.

To create another dataset for testing, a similar subsampling procedure was repeated for each volunteer. For testing data, the patches were not randomly subsampled but regularly in 5-pixel strides. As a result, all the pixels were subsampled at least once and most of the pixels (including the pixels of the brain area) were subsampled six or seven times.

CNN training and generating GenMVI volumes

The designed CNN network consists of a structure image block for acquiring local information from R1, R2, and PD maps, and a reconstruction block for adding that information to the corresponding SyMVF map. Details are described in Figure 1.

Fig. 1

The network architecture designed for this study. The network consisted of a structure image block and a reconstruction block. The structure image block has a contracting pathway on the left side and an expanding pathway on the right side. The block was designed to extract local information from the R1, R2, and PD maps. The reconstruction block reconstructs the corresponding synthetic myelin volume fraction (SyMVF) map into a new map image (generated myelin volume index: GenMVI). The numbers indicated at each convolutional (Conv) or transpose-convolutional (Conv-tr) layer describe the size of the input and output images (rows, columns, and channels) for the layer.

To create a GenMVI map for a particular volunteer, data from the other 19 were used to train a CNN for that volunteer (leave-one-out cross-validation; to lessen the effect of the small sample size). The 19 volunteers were randomly assigned to two groups of 15 and four subjects to create training and validation data sets. Gaussian noise was added to the training data to avoid overfitting.

Training was performed using thethe Tensorflow-GPU (Google, Mountainview, USA) (version 1.8.0)[28] platform with Keras[29] (version 2.1.6). Calculations were performed by computer equipped with: dual CPU, Xeon (Intel, Santa Clara, USA) E5-2623v4; dual GPU, Titan X Pascal (Nvidia, Santa Clara, USA), 12 GB GDDR5X; 128 GB random access memory; Ubuntu 16.04LTS. R1, R2, and PD patches were input to the structure image block, and the SyMVF patches were input to the reconstruction block. As loss-function, the root-mean-square-error (RMSE) with respect to MTMVI was obtained for both main output and auxiliary outputs, and then summed after multiplying RMSE of the auxiliary output by 0.2. The Adam algorithm[30] was applied for optimization, where the learning rate started from 0.0001 at the first epoch and then decreased according to the hyperbolic function, where lr(n) is the learning rate for epoch number n. The maximum number of epochs for training was set at 10, but training was aborted when the loss obtained at the end of each epoch using the validation data did not decrease for three consecutive epochs.

After completion of training, a GenMVI map was generated for each volunteer by test dataset. The dataset for each volunteer was input to a trained CNN (i.e. trained by the other 19 volunteers for each volunteer) to obtain main output as GenMVI patches. The output patches were re-orientated to form a whole GenMVI volume. For the pixels included in two or more patches, the median of the estimated values was applied.

Statistical analysis

Analysis based on averaged value of pixels inside ROIs

The 164 local ROIs were separately applied to the MTMVI, SyMVF, and GenMVI maps of the volunteers. The values of the pixels included in each ROI were averaged and recorded. The absolute difference of the averaged values of the SyMVF and MTMVI maps were calculated for each local ROI (ΔSy), and the same was done for the GenMVI and MTMVI maps (ΔGen). ΔSy and ΔGen were statistically compared for each of the following four regions: (i) cortical GM (consisting of 108 ROIs), (ii) subcortical GM (consisting of eight ROIs), (iii) WM (consisting of 48 ROIs), and (iv) whole brain (consisting of all 164 ROIs). The Wilcoxon signed-rank test was applied and P < 0.05 was considered significant.

For further comparison, Pearson’s correlation analysis was performed between MTMVI and SyMVF, and between MTMVI and GenMVI for the averaged values calculated for a total of 3280 local ROIs (164 ROIs from each of the 20 volunteers).

Pixel-based comparison within ROIs

ROIcGM, ROIsGM, ROIWM, and ROIWB were applied to each volunteer. The pixels in these four ROIs were extracted and used to calculate a pixel-based Pearson’s correlation coefficient for both SyMVF and GenMVI in comparison to MTMVI, for each volunteer. Distributions of SyMVF-based and GenMVI-based correlation coefficients for the 20 volunteers were compared statistically for all four ROIs. Wilcoxon signed-rank test was applied for this purpose and P < 0.05 was considered significant.

In addition, a similar pixel-based comparison was performed for ROICC.

Results

Examples of MTMVI, SyMVF, and GenMVI maps from the same volunteer are shown in Figure 2.

Fig. 2

Typical magnetization transfer-based myelin volume index (MTMVI), synthetic myelin volume fraction (SyMVF), and generated myelin volume index (GenMVI) maps obtained from the same volunteer. Visually, the contrast of the GenMVI map is more similar to MTMVI than that of the SyMVF map (See the corpus callosum area as a particular example).

Results based on averaged value of pixels inside ROIs

The distribution and relationship between the averaged values of the MTMVI, SyMVF, and GenMVI maps for the 164 local ROIs are illustrated in scatterplots (Figure 3). The median and minimum to maximum ranges of these metrics for cortical GM, subcortical GM, WM, and whole brain areas are indicated in Table 1. The median and range of the absolute differences with respect to MTMVI (i.e. ΔSy and ΔGen) are also shown in Table 1. The median value of ΔGen was smaller than that of ΔSy for all areas, with all differences being significant (P < 0.001). Pearson’s correlation coefficient obtained from all local ROIs from the 20 volunteers was larger for MTMVI and GenMVI (R = 0.86) than for MTMVI and SyMVF (R = 0.77) (Figure 3). Both correlations were significant (P < 0.001).

Fig. 3

Scatter plots of the average values in the 164 local ROIs: correlation between (left) the magnetization transfer-based myelin volume index (MTMVI) and synthetic myelin volume fraction (SyMVF) maps (R = 0.77), and (right) the MTMVI and generated myelin volume index (GenMVI) maps (R = 0.86). Overall, the values are more consistent between maps for the latter comparison.

Table 1

The absolute errors between the SyMVF and MTMVI maps (ΔSy), and between the GenMVI and MTMVI maps (ΔGen)

	Cortical GM	Subcortical GM	WM	Whole brain
MTMVI	0.23 [0.08, 0.39]	0.24 [0.19, 0.30]	0.32 [0.17, 0.49]	0.25 [0.08, 0.49]
SyMVF	0.10 [0.06, 0.15]	0.20 [0.14, 0.28]	0.32 [0.19, 0.43]	0.11 [0.06, 0.43]
GenMVI	0.24 [0.14, 0.35]	0.25 [0.21, 0.31]	0.32 [0.20, 0.45]	0.26 [0.14, 0.45]
ΔSy* 1	0.13 [0.00, 0.29]	0.04 [0.00, 0.11]	0.02 [0.00, 0.12]	0.10 [0.00, 0.29]
ΔGen* 1	0.02 [0.00, 0.15]	0.01 [0.00, 0.06]	0.01 [0.00, 0.11]	0.02 [0.00, 0.15]
P-value* 2	<0.001	<0.001	<0.001	<0.001

Comparison between ΔSy and ΔGen (Wilcoxon signed-rank test, P < 0.05 considered as significant). Results are presented as median [min, max]. GM, gray matter; WM, white matter; MTMVI, myelin volume index calculated from magnetization-transfer images; SyMVF, myelin volume fraction estimated with SyMRI (name of commercial software); GenMVI, generated myelin volume index; ΔSy, |SyMVF − MTMVI|; ΔGen, |GenMVF − MTMVI|.

Results of pixel-based comparison within ROIs

Distributions of pixel-based correlation-coefficients for the 20 volunteers for each of ROIcGM, ROIsGM, ROIWM, and ROIWB are illustrated in box-plot graphs (Figure 4). Strong correlation between MTMVI and GenMVI was found for all four ROIs, where median values of the correlation coefficients were always higher than 0.80. On the other hand, those for SyMVF were lower than 0.7 except for ROIWM, for which the value was 0.70. Distributions for SyMVF and GenMVI were significantly different for all four ROIs (P < 0.001).

Fig. 4

The box plots show the distribution of pixel-wise correlation coefficients obtained by comparing the synthetic myelin volume fraction (SyMVF) and generated myelin volume index (GenMVI) maps with the magnetization transfer-based myelin volume index (MTMVI) map for the four ROIs corresponding to the cortical GM, subcortical GM, WM, and whole brain (i.e. ROIcGM, ROIsGM, ROIWM, and ROIWB), for all 20 volunteers. The median values are higher for the GenMVI map than for the SyMVF map, and the differences of the distributions are significant in all four areas (Wilcoxon signed-rank test, P < 0.001).

The result of an additional pixel-wise comparison for ROICC is also illustrated in a box-plot graph (Figure 5). The correlation was moderate for GenMVI (median value 0.56), but it was stronger than that for SyMVF (median value 0.21). The difference in distributions was significant (P < 0.001).

Fig. 5

The box plots show the distribution of pixel-wise correlation coefficients obtained by comparing the synthetic myelin volume fraction (SyMVF) and generated myelin volume index (GenMVI) maps with respect to the magnetization transfer-based myelin volume index (MTMVI) map using the corpus-callosum ROI (ROICC). The correlation coefficient obtained for a GenMVI map is always higher than that obtained for the corresponding SyMVF map. The difference is significant (Wilcoxon signed-rank test, P < 0.001).

Discussion

This study was designed to evaluate the usefulness of a CNN-based method for estimating a MVI map from RSRI images. Combining the overall results, it seems that the method proposed in this study successfully reconstructs the contrast of a SyMVF map of healthy brain into a new contrast that is more strongly related to the corresponding MTMVI map.

CNN architecture

The structure image block aimed to compensate for the potential weak point of a SyMVF map, namely, not containing any information about local tissue structure. The structure image block consists of many convolutional layers between input and output, which means that the value assigned to a pixel at the output incorporates information from a relatively wide area surrounding the pixel in the input image (i.e. maximum 32 × 32-pixel area). This block was designed based on U-net,[16] which is a network that has achieved great success when used to segment various anatomical and histological images.[31-33]

The overall CNN designed in this research aimed to utilize the high capability of CNN for flexible image segmentation and reconstruction, while at the same time designing the network so that the priority of the SyMVF map as input is relatively high (e.g. shortcuts in reconstruction block that sent SyMVF images to later layers without being merged with the input from the structure image block). The design aimed to lessen the black box problem[17,34] as much as possible by using SyMVF, which has a logical foundation (i.e. Bloch simulation), as basic starting point and giving more priority in the overall network design.

Comparison of SyMVF and GenMVI with respect to MTMVI

Visually, the contrast of the GenMVI map was closer than the contrast of SyMVF to that of MTMVI (Figure 2). In particular, characteristics of the corpus callosum were better reproduced in the GenMVI map than in the SyMVF map. The corpus callosum was where the contrast was especially different between the SyMVF and MTMVI maps in a previous study.[3] The method proposed in this study appears to have improved this problem.

From analysis using the averaged values of the 164 local ROIs, values of the GenMVI map were more similar to the MTMVI map for different brain areas than those of the SyMVF map (Figure 3 and Table 1). In addition, even though the correlation for SyMVF was also strong (R = 0.77), which is consistent with a previous study,[3] the overall linear correlation with MTMVI was higher for GenMVI (R = 0.86) (Figure 3).

The results of the pixel-based comparisons further support the results of the atlas-based comparison, namely, that the GenMVI map has a stronger correlation with MTMVI than the SyMVF map does (Figure 4). Pixel-based comparison was added because the results using averaged values for each ROI are potentially biased by a possible difference in pixel number in each ROI.

The pixel-based comparison performed for the corpus callosum indicated higher correlation for GenMVF over SyMVF, with a statistically significant difference between the two distributions (P < 0.001) (Figure 5). This result is consistent with the visual evaluation of the region (Figure 2).

Currently, one of the great advantages of SyMVF over simpler and less time-consuming myelin-related indices may be its high correlation with the widely used MTMVI map.[3] The proposed GenMVI maps may improve on the advantages of the SyMVF maps to expand the capability toward possible clinical use. However, MTMVI is still not a golden standard that is comparable to histology. Thus, the fact that the contrast of the GenMVI map was closer than the SyMVF map to the MTMVI map does not directly mean that the GenMVI map is more accurate than the SyMVF map. Further study using histology-based measures of myelin volume as a target during training is desirable. Moreover, since the GenMVI map has the limitation that its logical foundation is relatively weak, further clinical validation is also important (See supplementary material for a preliminary clinical application. Supplementary material is available online).

A possible improvement for the proposed method would be to include new information in addition to the current input. The most promising candidates are various diffusion-related metrics,[35,36] such as apparent diffusion coefficient, fractional anisotropy (FA), and other metrics obtained from diffusional kurtosis imaging. It is well established that these parameters are closely related to the local micro-structure and myelin volume content in some tissues including demyelinating lesions.[35] The additional information may help to appropriately estimate myelin volume especially when the target will be expanded to pathological brains.

There are some limitations in this study. First, the MTMVI maps were warped to register the images to R1, R2, PD and SyMVF maps obtained from RSRI. A small mis-registration might have affected the training, as well as the final results of the statistical analysis. Second, the possibility of overfitting has not been completely ruled out, even though we used cross-validation to be strict to this problem as much as possible. This is because the subjects were all healthy volunteers, and the dispersion of their myelin measure were relatively small. This also means that it is still not evident that this estimation is adequate for the patients (e.g. patients with some white matter disease) as similarly for the healthy volunteers. Including patients and increasing the number of the subject in addition may help overcoming this limitation. Nevertheless, clinical validation is necessary as described above.

Third, The CNN architecture of this study is not proven as the best for the task. The hyperparameters of the CNN might remain some room for further optimization. In addition, the structure image block, which was designed based on U-net in this study, might become more effective by changing the main structure to other established ones such as SegNet.[37] However, we did not explore further in this theme because such optimization might better be done with large sample size.

Conclusion

In conclusion, the deep-learning-based method proposed in this study generated a MVI from RSRI that incorporates more specific information about local tissue properties than the existing technique. This method might be useful for the future clinical applications; however, further work is necessary to validate the proposed method because this study was performed only for healthy volunteers.