Simultaneous Mapping of Water Diffusion Coefficients and Metabolite Distributions of the Brain Using MR Spectroscopic Imaging Without Water Suppression

Objective: To simultaneously map water diffusion coefficients and metabolite distributions of the brain in magnetic resonance spectroscopic imaging (MRSI) experiments within a clinically feasible time. Methods: A diffusion-preparation module was introduced in water-unsuppressed MRSI acquisition sequence to generate diffusion weighting of the water signals. Fast spatiospectral encodings were achieved using echo-planar spectroscopic imaging readouts with blipped phase encodings for sparse sampling. Navigator signals were embedded in the data acquisition sequence, which were used for detection of data corrupted by physiological motion in the diffusion preparation period. In data processing, a novel model-based method was developed to effectively use sparse (k, t)-space spectroscopic signals for reconstruction of the spatial distributions of water diffusion coefficients and metabolite concentrations. Results: Both phantom experiments and in vivo experiments were carried out to evaluate the feasibility and performance of the proposed method. In an 8-minute scan, diffusion weighted images and apparent diffusion coefficients map at 2.0×1.0×1.0 mm3 were obtained simultaneously with metabolite maps at 2.0×3.0×3.0 mm3 nominal resolution. Conclusion: We demonstrated the feasibility of using the unsuppressed water signals from MRSI experiments to map the water diffusion coefficients of brain tissues and proposed a novel method to achieve simultaneous mapping of water diffusion coefficients and metabolite distributions. Significance: The proposed method provides a unique imaging tool for simultaneous diffusion and metabolic imaging. This method is expected to be useful for various brain imaging applications.


I. INTRODUCTION
M R DIFFUSION imaging and spectroscopic imaging (SI) are two key imaging tools for brain mapping. Diffusion MRI is often used to map tissue microstructural properties [1], [2], [3], [4] while MRSI provides spatial maps of brain metabolites and neurotransmitters [5], [6], [7], [8]. Recently, there is an increasing interest in including both diffusion MRI and MRSI in disease studies [9], [10], [11], [12]. It has been shown that combining information from both imaging modalities can improve the prediction of stroke outcome [13], increase the accuracy of cancer detection [14], and improve the differentiation between benign and malignant tumors [15].
In the existing studies, diffusion MRI and MRSI data were acquired separately using different pulse sequences. More specifically, pulsed-gradient spin-echo echo-planar imaging (EPI) has been widely used for diffusion MRI due to its high acquisition speed and relative immunity to physiological motion [16]. However, the EPI trajectory is sensitive to B 0 field inhomogeneity and T 2 * decay, which can lead to significant geometric distortions and signal loss [2]. Phase encoding-based MRSI sequences have been used as standard MRSI methods for decades. But due to the low concentration of metabolites in the brain, these traditional MRSI methods have long been limited by low signal-to-noise ratio (SNR), low resolution on the order of centimeters, and long scan time more than 20 minutes [5], [17]. Acquiring diffusion MRI and MRSI separately using conventional sequences often leads to a long total scan time. Another important technical issue lies in the challenge of accurately registering images from two acquisition modalities with different image characteristics (e.g., geometric distortions often associated with diffusion MRI data and the low SNR and low-resolution nature of MRSI data).
In this work, we present a new method for simultaneous diffusion and metabolic imaging by leveraging these advances. The proposed method keeps the basic features of the SPICE sequence for acquisition of water-unsuppressed MRSI data and introduces a diffusion-preparation module to generate diffusion weighting of water spectroscopic signals. The (k, t)-space of diffusion-weighted signals is covered very sparsely using CAIPIRINHA (Controlled Aliasing in Parallel Imaging Results in Higher Acceleration) trajectories [48]. In data processing, a generalized-series (GS) model incorporating the unsuppressed water signals is used for high-quality reconstruction and correction of the data corrupted by physiological motion [49]. As a result, distortion-free diffusion weighted images and apparent diffusion coefficient (ADC) maps at 2.0×1.0×1.0 mm 3 as well as high-quality metabolite maps at 2.0×3.0×3.0 mm 3 nominal resolution were simultaneously obtained in an 8-minute scan. Phantom experiments were performed on the NIST phantom for validation, and in vivo experiments were carried out on both healthy subjects and tumor patients to demonstrate the performance and potential of the proposed method.

A. Data Acquisition
The proposed sequence keeps the basic features of the SPICE sequence for acquisition of water-unsuppressed MRSI data and includes a diffusion-preparation module for mapping of water diffusion, as illustrated in Fig. 1(a). As the basic SPICE method, the proposed sequence eliminates water suppression pulses, acquires free induction decay (FID) signals in echo planar spectroscopic imaging (EPSI) trajectories, and uses ultra-short TE (echo time, 1.6 ms) and short TR (repetition time, 160 ms) to achieve good SNR and high scanning speed. These features enable fast high-resolution metabolite mapping and simultaneously generate water spatiotemporal signals without diffusion weighting. To generate diffusion weighting of the water signals, diffusion preparation pulses are added before excitation. The diffusion preparation module contains a pair of nonselective 90-degree pulses and two pairs of nonselective 180-degree refocusing pulses, with the diffusion gradients being placed between these RF pulses. This diffusion preparation scheme is designed to reduce the sensitivity to B 0 and B 1 inhomogeneities, as in the previous method used for T 2 -preparation [50]. In addition, three The basic SPICE sequence was extended with a diffusion-preparation module for simultaneous diffusion and metabolic imaging. The diffusionpreparation module included three pairs of hard pulses for insensitivity to field variation. At the end of each TR, linear navigators were embedded to collect signals for correction of physiological motion induced error. (b) Sampling pattern of the proposed method in (k, t)-space. Instead of being fully sampled, (k, t)-space of both SPICE frame and diffusion frames were sparsely sampled using spatiotemporal CAIPIRINHA trajectories.
linear navigators in orthogonal directions are embedded at the end of each TR to detect whether the data in the TR is corrupted by physiological motion. However, the EPSI trajectories in the basic SPICE sequence have relatively low spatial encoding efficiency; using them to sample the (k, t)-space of multiple frames with different diffusion weightings would lead to a long scan time (e.g., around 30 minutes for each frame). Therefore, in the proposed method, (k, t)-space is sparsely sampled in variable density to achieve high imaging speed, as illustrated in Fig. 1(b). More specifically, the metabolite signals and high-resolution water signals without diffusion weighting are acquired as in the basic SPICE sequence (denoted as SPICE frame), which is followed by the acquisition of frames with different b-values and directions (denoted as diffusion frames). The metabolite signals are fully sampled in central k-space to preserve SNR, while the water signals in both SPICE frame and diffusion frames are highly sparsely sampled in spatiotemporal CAIPIRINHA trajectories.
In the current implementation of our proposed method, the sequence parameters are set as follows: FOV = 240×240×72 mm 3 ; TE = 1.6 ms; TR = 160 ms; echo space = 1.76 ms; flip angles = 27˚/60˚(for SPICE frame/diffusion frames); matrix size of water signals: 120×216×72 (corresponding to 2.0×1.0×1.0 mm 3 resolution); matrix size of metabolite signals: 120×78×24 (corresponding to 2.0×3.0×3.0 mm 3 resolution). With this setup, acquisition of the SPICE frame takes 7 minutes [46]. In covering (k, t)-space of the diffusion frames, the ky-dimension is uniformly under-sampled by a factor of 3 (except for the fully sampled central 10 ky lines), while the (kz, t)-space is sparsely undersampled by a factor of 72 using CAIPIRINHA trajectories, resulting in a 10-second scan time for each diffusion frame. Therefore, a typical simultaneous MRSI and diffusion MRI scan including one SPICE frame and four diffusion frames (one frame with zero b-value and three frames with orthogonal diffusion directions) can be scanned within 8 minutes.

B. Image Reconstruction
The MRSI data in the SPICE frame were processed using the original SPICE pipeline to generate metabolite maps and highresolution water signals without diffusion weighting as previously reported [40], [41], [42], [43], [44], [51]. The main tasks to process the diffusion imaging data include: (1) reconstruction of the spatiospectral functions from the sparse measurements, and (2) correction of the effects induced by physiological motion.
To solve the sparse reconstruction problem, a GS model incorporating the high-resolution unsuppressed water signals from the SPICE frame was used. We took advantage of the fact that the signals collected in the diffusion frames and the SPICE frame shared the same acquisition sequence except for the diffusion-preparation module. In other words, we treated them to have same temporal evolution process except for the diffusion weighting. Therefore, the relationship between the spatiotemporal functions of the water signals in the SPICE frame (denoted as ρ r (x, t)) and those in the diffusion frames (denoted as ρ i (x, t)) could be expressed using a low order GS model [49]: where i denotes a diffusion frame with a specific b-value and direction, c n (x, i) the GS model coefficients, N i the GS model order, and Δf a frequency step. Then reconstruction of ρ i (x, t) from the sparsely measured data (denoted as d i ) became solving the following optimization problem frame by frame: where c i and ρ r are the vector forms of {c n (x, i)} and ρ r (x, t), Ω, F, B, S, G, W represent operators associated with k-space sampling, Fourier transform, field inhomogeneity effect, coil sensitivity weighting, GS modeling and edge-weighted total variation, respectively. λ is the regularization parameter, and this optimization problem was solved using the conjugate gradient method. After the GS coefficientsĉ i were determined, the reconstructed signals of each diffusion frameρ i (x, t) were synthesized using the GS model and the reference SPICE water signals, as shown in (1).

C. Motion Correction
With the tip-up pulse in the diffusion-preparation module, the phase error induced by physiological motion can cause errors in both magnitude and phase [52]. To correct these errors, we first identified the corrupted TRs and their corresponding (k, t)-space data using the navigator signals. Then, the corrupted data were replaced by synthetic data generated using the uncorrupted data by a model-based interpolation. The final reconstruction results were generated using the corrected data.
More specifically, the linear navigators in three directions collected three orthogonal k-space lines passing through the k-space center. These three navigator signals acquired in each TR were combined into one signal via sum-of-square. Since physiological motion caused magnitude loss (so the navigators in the corrupted TRs had lower signal amplitude), the navigator with the highest overall magnitude was used as the reference navigator of the uncorrupted TRs. Then the navigators of the corrupted TRs were identified as those with low correlation with the reference navigator. After identifying the corrupted TRs and the corresponding (k, t)-space data, the remaining uncorrupted (k, t)-space data were used to regenerate and replace the corrupted data using interpolation. The interpolation could be achieved by frame-by-frame reconstruction as shown in (2), only using the uncorrupted (k, t)-space data. But when the amount of motion-corrupted data is large, the amount of data used for reconstruction would become very small, making the reconstruction problem ill-conditioned. Thus, an additional low-rank constraint across all the frames was imposed to include more data from multiple frames for reconstruction. More specifically, the GS coefficients of all the diffusion frames {c n (x, i)} were expressed using a low-rank subspace model: where {φ n,l (i)} is the basis function across the diffusion frames and {b n,l (x)} the corresponding spatial coefficients. The basis function was estimated from one set of training data acquired using the proposed sequence, but with fully sampled (k, t)space in low resolution. More specifically, the matrix size of the sequence for training data was 120×78×24 (corresponding to 2.0×3.0×3.0 mm 3 resolution, fully sampled), while other parameters were kept the same. With this setup, each frame took 5-minute scan time and 5 frames (including one SPICE frame, and four diffusion frames with b-value = 0, and 1000 s/mm 2 in three orthogonal directions) were acquired on each subject. In data processing, first, the GS coefficients were obtained from the training data by fitting to (1) frame by frame. Then, the basis function was estimated from the GS coefficients using singular value decomposition. With the basis function predetermined, the reconstruction using the uncorrupted data of all the diffusion frames became: where b and Φ are tensor forms of {b n,l (x)} and {φ n,l (i)}, d u and Ω u represent the vector form of the uncorrupted (k, t)space data and the corresponding (k, t)-space sampling mask, respectively. Similar to (2), λ is the regularization parameter, and the optimization problem was also solved using the conjugate gradient method. In both equations, the λ was chosen as very small values to impose only weak spatial constraint. After the reconstruction, synthetic (k, t)-space data were used to replace the corrupted data. Then this set of corrected data was used to reconstruct the final diffusion images using (2).

D. Experimental Setup
The proposed pulse sequence was implemented on 3T Siemens Prisma/Skyra scanners (Siemens Medical Solutions, Erlangen, Germany). To demonstrate the feasibility and potential of the proposed method, both phantom and in vivo experiments were carried out. In the phantom scans, the NIST/ISMRM diffusion phantom was used. To validate the reconstruction method, a set of full data was acquired using the proposed sequence without sparse sampling (b-value = 0 s/mm 2 , and 500 s/mm 2 in three orthogonal directions; scan time = 42×5 min); then the sparse data were simulated by retrospective sampling. Besides, a typical EPI diffusion sequence (FOV = 240×240×72 mm 3 ; TE = 128 ms; TR = 6700 ms; matrix size = 160×160×36; resolution = 1.5×1.5×2.0 mm 3 ; echo space = 1.01 ms, b-value = 0 s/mm 2 , and 500 s/mm 2 in three orthogonal directions, scan time = 1:30 min) was also scanned for comparison. In vivo studies were performed on healthy subjects and tumor patients, which were all approved by the local Institutional Review Board. A written consent form was obtained prior to each study. The scan protocol included a typical MPRAGE imaging (FOV = 240×240×192 mm 3 , TE = 2.29 ms, TI = 900 ms, TR = 1900 ms, matrix size = 256×256×192 (resolution = 1.0×1.0×1.0 mm 3 ), scan time = 5 min), simultaneous diffusion imaging and MRSI using the proposed method (b-value = 0, 500 or 1000 s/mm 2 in three orthogonal directions; scan time = 8 min), and a diffusion imaging scan using the EPI diffusion sequence (FOV = 240×240×72 mm 3 ; TE = 128 ms; TR = 6700 ms; matrix size = 160×160×36; resolution = 1.5×1.5×2.0 mm 3 ; echo space = 1.01 ms, b-value = 0, and 1000 s/mm 2 in three orthogonal directions, scan time = 1:30 min) for comparison. For the tumor patients, contrast agent (gadolinium diethylenetriaminepentaacetic acid) was intravenously administered for contrast enhancement of the MPRAGE images as part of the clinical protocol. Fig. 2 shows a set of diffusion weighted images and ADC maps of the NIST diffusion phantom obtained by the proposed method and the traditional EPI diffusion method. Compared with the reference image (water image obtained in the SPICE frame), the shape of several phantom vials was distorted in the EPI images due to susceptibility effects. The diffusion images of the proposed method showed little geometric distortion and the phantom vials were in good circular shape. Note also that the results reconstructed from the sparse data (using retrospective sampling) showed good agreement with those from the full data, which indicated the effectiveness of the reconstruction method. Fig. 3(a) shows the comparison of ADC values measured by the proposed method and the typical EPI diffusion sequence, and Fig. 3(b) shows the comparison between the fully sampled  data and the retrospectively sampled sparse data. In these comparisons, mean ADC values of the phantom tubes measured using different methods were calculated and displayed as the scatter plots. The regression lines of the scatter points were also computed and compared with the identical line, and the Pearson's correlation coefficients were calculated. As shown in Fig. 3, the ADC values obtained by the proposed method had good consistency with those obtained using the typical EPI diffusion sequence (γ = 0.9851). The results obtained from the sparse data were in good agreement with those from the fully sampled data (γ = 0.9998), which indicated the reconstruction method did not introduce noticeable bias in the ADC values. Fig. 4 shows a set of diffusion images acquired from a healthy subject, including diffusion-weighted images of different bvalues, diffusion-weighted images of different TEs and derived ADC maps. From Fig. 4(a), we can see high-resolution diffusion weighted images of the brain were successfully reconstructed from the highly sparse data. Since the diffusion images and the water images of the SPICE frame were acquired from the same EPSI trajectories, they were spatially well aligned. Each diffusion frame contains 148 images of different TEs, which are of different contrasts, as some of them shown in Fig. 4(b). The  ADC map obtained using the traditional EPI method was also shown in Fig. 4(c) for comparison. As can be observed, there were noticeable geometric distortions near the frontal sinus regions in the ADC map produced by the EPI method, which were not present in the images obtained using the proposed method. Comparing the ADC values obtained using the proposed method and the EPI method, a good agreement was found. The mean ADC values of gray matter and white matter were 0.959×10 −3 mm 2 /s, 0.734×10 −3 mm 2 /s in the EPI data and 0.941×10 −3 mm 2 /s, 0.752×10 −3 mm 2 /s using the proposed method. Fig. 5 shows a representative set of results including diffusion weighted image (b = 1000 s/mm 2 ) and ADC map at 2.0×1.0×1.0 mm 3 , and metabolite maps (including NAA, Cr, and Cho maps) at 2.0×3.0×3.0 mm 3 , obtained using the proposed method in a total 8-minute scan. All these images were generated simultaneously from the same acquisition, thus were naturally co-registered. The right panel shows several spatially localized spectra, illustrating the high quality of the reconstructed spatiospectral functions. Another set of results acquired from a patient diagnosed with high-grade glioma (WHO IV) are shown in Fig. 6. The position of tumor could be observed on the contrast-enhanced T 1 -weighted image. Compared with the normal tissues, hyperintensity in the diffusion images, increased ADC values, decreased NAA and increased Cho could be observed in the tumor, which were consistent with the literature reports [53]. The alteration in tissue diffusion properties, loss of neurons, and abnormal membrane turnover in tumor were all reflected from the imaging results obtained by the proposed method in an 8-minute scan.

IV. DISCUSSION
This paper presents a novel method for simultaneous diffusion imaging and metabolic imaging using SPICE. The SPICE sequence was extended with a diffusion-preparation module for 3D distortion-free diffusion imaging. Besides, highly sparse sampling was employed for acceleration, making the total scan time clinically feasible.
The SPICE sequence used for MRSI acquisition in the proposed method has several unique features including FID acquisition, no water suppression, rapid EPSI trajectories, ultrashort TE, short TR, and sparse sampling to enable fast high-resolution metabolite mapping, providing significantly improved efficiency compared with traditional MRSI methods [39], [43], [46]. Besides, the unsuppressed water signals not only facilitate correction of system imperfections, but also enable simultaneous mapping of other imaging contrast, such as functional imaging [45], susceptibility mapping [46], [47], and diffusion imaging proposed in this work.
Compared with the widely used EPI diffusion sequence, the proposed sequence has several advantages for diffusion imaging: 3D volumetric excitation and acquisition (thus, higher SNR and better slice profile [2]), free of geometric distortion associated with EPI trajectories, and flexibility for multi-shot acquisition with potential for higher-resolution diffusion imaging. In sampling (k, t)-space, typical EPI trajectories encode different ky lines at different t and combine them to form an image, which yields high spatial encoding efficiency but leads to geometric distortion due to B 0 field inhomogeneity and T 2 * decay effects [2], [3]. EPSI trajectories used in the proposed method cover ky and kz dimensions using phase encoding, which means full 3D data at each t in (k, t)-space are acquired, thus significantly reducing the geometric distortion associated with EPI trajectories. EPSI trajectories has relatively low spatial encoding efficiency for diffusion imaging, while the associated acceleration issue was well addressed in our method using sparse sampling and a GS model-based reconstruction. The GS model used in proposed method incorporates the water spatiotemporal functions of the SPICE frame as a reference and transforms the reconstruction issue into estimation of GS model coefficients. This GS model significantly reduces the degrees of freedom so makes the reconstruction from such highly sparse data possible.
The currently implemented diffusion-preparation module includes three pairs of hard pulses to reduce the sensitivity to B 0 and B 1 inhomogeneities. One alternative is to use adiabatic pulses which are insensitive to field variations, but they have issues of high energy deposition, especially when the TR is as short as 160 ms. Fig. 7 compares the currently used preparation pulses with typical one-refocusing scheme using hard pulses and twice-refocusing scheme using two adiabatic pulses. As can be seen, the proposed preparation pulses have better insensitivity to field variations than one-refocusing hard-pulse scheme and similar performance to the twice-refocusing adiabatic-pulse scheme. It should be also noted that the proposed diffusion-preparation module contains four refocusing pulses and multiple diffusion gradients with alternating polarities, which can reduce the cumulated eddy current effects caused by the diffusion gradients [54], [55], [56]. In many other diffusion-preparation methods, a crusher gradient called stabilizer gradient was used to address the error induced by physiological motion but at the price of half the SNR [52], [57], [58]. No stabilizer gradient was used in the proposed method. The proposed method utilized the SPICE frame, GS modeling, and subspace modeling to interpolate and replace the corrupted (k, t)-space data, which corrected the motion-induced errors effectively, thus the stabilizer was discarded for better SNR. Fig. 8 shows the effectiveness of the proposed method in correcting motion-induced errors. Without correction, there were noticeable artifacts in the reconstructed images. By excluding the corrupted TRs, the artifacts were significantly reduced but remained noticeable; after correction using the proposed subspace-model based method (as in (3) and (4)), the artifacts became negligible and image features were well preserved. However, the current correction method could only work under the assumption that the corrupted TRs were only a small portion (<20%) of the total measurements, which was the case in most of our experiments.
Several potential improvements can be made in future development. First, the current diffusion preparation pulses were relatively insensitive to field variations but there was still signal loss when large field inhomogeneity was present. Utilization of SAR-efficient adiabatic pulses may further improve the field insensitivity. Second, current method to correct the motion induced error may not work well when a large portion of the measurements were corrupted, which would limit its applications in higher b-values or more diffusion directions, like diffusion tensor imaging. More effective correction methods with 2D/3D navigators may further improve the robustness of the proposed method [25], [59]. Third, current choice of regularization is an L2-norm based regularization. It can be replaced with L1 or nuclear norm to provide stronger spatial constraint, but at the expense of computational efficiency. The data collected in the proposed method have high dimensionality, including three spatial dimensions, one temporal/spectral dimension, one coil/channel dimension, and one diffusion dimension. Even in the linear formulation with L2 regularization, the computational burden is still quite high. Therefore, more effective methods to rapidly solve tensor problem will be very helpful. Moreover, most of the operations in our method are formulated into matrix operations, which do not fully utilize the power of the low rankness of the tensor structure. More effective methods to solve the tensor problem could further improve our method.

V. CONCLUSION
This paper presents a new method for simultaneous diffusion imaging and metabolic imaging of the brain, which is able to achieve distortion-free diffusion imaging at 2.0×1.0×1.0 mm 3 resolution and metabolic imaging at 2.0×3.0×3.0 mm 3 nominal resolution. With further development, the proposed method may provide a useful tool for various clinical applications.