A fundamental limitation with the conventional one-dimensional (1D) proton MRSI (1H-MRSI) is that most of the detectable resonances lie in a small chemical shift range of approximately 4 ppm. As a result, the spectra of these metabolites often overlap with each other (e.g., Glu and GABA), making their accurate detection difficult. This limitation can be overcome by using multidimensional spectroscopy, e.g., two-dimensional (2D) J-resolved spectroscopy, where resonances are distributed in multidimensional spaces with increased spectral dispersion. However, the additional dimension of spectral information is obtained at the cost of increased data acquisition time, limiting its practical utility. Fast imaging sequences and compressed sensing based methods have been proposed to accelerate multidimensional 1H-MRSI.^1-4 This work presents a novel tensor-based approach to accelerated high-resolution multidimensional 1H-MRSI.

To take advantage of data correlation in multiple dimensions for accelerated imaging, we represent multidimensional 1H-MRSI signals as high-order partially separable (HOPS) functions:⁵ $$\rho(\mathbf{x},t_1,t_2)=\sum_{l=1}^{L} \sum_{m=1}^{M} \sum_{n=1}^{N} c_{lmn} \theta_{l}(\mathbf{x}) \phi_{m}(t_1) \psi_{n}(t_2),$$where $$$\{\theta_{l}(\mathbf{x})\}_{l=1}^{L}$$$, $$$\{\phi_{m}(t_1)\}_{m=1}^{M}$$$, and $$$\{\psi_{n}(t_2)\}_{n=1}^{N}$$$ are basis functions that describe the variation of the spatiotemporal function $$$\rho(\mathbf{x},t_1,t_2)$$$ along the spatial ($$$\mathbf{x}$$$), indirect evolution (due to J-coupling, $$$t_1$$$), and direct evolution (due to chemical shift, $$$t_2$$$) axes, respectively, and $$$\{ c_{lmn} \}_{l,m,n=1}^{L,M,N}$$$ are the corresponding coefficients. Mathematically, the HOPS model can be viewed as low-rank tensors. In fact, after discretization, $$$\rho(\mathbf{x},t_1,t_2)$$$ can be represented by an order-3 tensor in the Tucker form with a tensor rank ($$$L$$$,$$$M$$$,$$$N$$$).⁶

The HOPS model enables the use of special data acquisition schemes to sparsely sample $$$(k,t_1,t_2)$$$-space. An example is shown in Fig. 1. Two “training” datasets ($$$D_{t_1}$$$ and $$$D_{t_2}$$$) are acquired with limited k-space coverage for determining $$$\{\phi_{m}(t_1)\}_{m=1}^{M}$$$ and $$$\{\psi_{n}(t_2)\}_{n=1}^{N}$$$, respectively. More specifically, $$$D_{t_1}$$$ samples $$$(k,t_1,t_2)$$$-space densely along $$$t_1$$$; and $$$D_{t_2}$$$ samples $$$(k,t_1,t_2)$$$-space densely along $$$t_2$$$ but sparsely along $$$t_1$$$. A separate sparsely sampled “imaging” dataset ($$$D_{2}$$$) with extend k-space coverage is acquired to recover the image function at high resolution.

In image reconstruction, we first estimate the subspace structure ($$$\{\phi_{m}(t_1)\}_{m=1}^{M}$$$ and $$$\{\psi_{n}(t_2)\}_{n=1}^{N}$$$) explicitly from the “training” datasets. To determine $$$\{\phi_{m}(t_1)\}_{m=1}^{M}$$$, a Casorati matrix is formed using the data in $$$D_{t_1}$$$, in which the $$$t_1$$$ coordinates of the measured data change with the rows. $$$\{\phi_{m}(t_1)\}_{m=1}^{M}$$$ are then determined by calculating the left singular vectors of the Casorati matrix using SVD. $$$\{\psi_{n}(t_2)\}_{n=1}^{N}$$$ are determined similarly. The effects of $$$B_0$$$ inhomogeneities are corrected before determining $$$\{\psi_{n}(t_2)\}_{n=1}^{N}$$$.^7,8

Finally, $$$\{ c_{lmn} \}_{l,m,n=1}^{L,M,N}$$$ and $$$\{\theta_{l}(\mathbf{x})\}_{l=1}^{L}$$$ are determined by fitting the data in $$$D_{2}$$$ using the estimated subspaces ($$$\{\hat{\phi}_{m}(t_1)\}_{m=1}^{M}$$$ and $$$\{\hat{\psi}_{n}(t_2)\}_{n=1}^{N}$$$)$$\min \parallel \mathbf{s}_{2}- \mathbf{F}_{B} \{ \sum_{l=1}^{L} \sum_{m=1}^{M} \sum_{n=1}^{N} c_{lmn} \theta_{l}(\mathbf{x}) \hat{\phi}_{m}(t_1) \hat{\psi}_{n}(t_2)\}\parallel_2^2 + R_{1}(\theta_{l}(\mathbf{x})) + R_{2}(c_{lmn}),$$ where $$$\mathbf{s}_{2}$$$ contains the sparsely sampled data in $$$D_2$$$, $$$\mathbf{F}_{B}$$$ is a Fourier encoding operator taking $$$B_0$$$ field inhomogeneity into account, the first regularization term is used to incorporate the prior knowledge of the spatial distributions of metabolites,⁸ and the second regularization term penalizes the sparsity of $$$ c_{lmn}$$$.

The proposed method has been validated by performing 2D J-resolved EPSI 1H-MRSI experiments on both phantoms and healthy subjects (approved by our local IRB) using a 3T Siemens Trio scanner.

Figures 2 and 3 show the results from equivalent-time acquisitions on a phantom with 8 vials filled with solutions of NAA, Cr, Cho, mI, Glu, Gln, and GABA. In the proposed method, the $$$D_{t_1}$$$ dataset was acquired with 6x12 spatial encodings, 20 TEs, and 256 echoes (0.78 ms echospacing); the $$$D_{t_2}$$$ dataset was acquired using the same sequence with 10x12 spatial encodings and 8 TEs; and the $$$D_2$$$ dataset was acquired with 64x64 spatial encodings, 20 TEs, and 128 echoes (1.98 ms echospacing). The $$$D_2$$$ dataset was randomly undersampled in the $$$(k,t_1)$$$ plane with a reduction factor of 4. The proposed method was compared with a low-resolution EPSI acquisition and the group sparsity method.³ As can been seen, the proposed method recovered the spatiospectral distribution of the phantom with both high resolution and SNR.

Figure 4 shows a set of in vivo results using the proposed method, demonstrating the promising capability of the proposed method in recovering a high-resolution, high-SNR multidimensional 1H-MRSI spatiospectral distribution. The data acquisition parameters were the same as in the phantom experiment. The nuisance signals were removed using the subspace-based method.⁹ The $$$D_2$$$ data were retrospectively undersampled after nuisance signal removal.

This paper presents a novel tensor-based method to achieve high-resolution multidimensional 1H-MRSI. The proposed method has been validated using phantom and in vivo experiments, producing encouraging results. It should enhance the practical utility of multidimensional 1H-MRSI.