Dynamic MRSI measures the temporal changes of metabolite concentrations by acquiring a time series of MRSI data. These data can be used in a range of application, including the study of the response of a metabolic system to a perturbation. A well-known example is the study of the depletion and recovery of PCr of the muscle before and after exercise using phosphorous (31P) NMR.¹ However, high-resolution dynamic MRSI is challenging due to poor SNR resulting from the low concentrations of metabolites. It is even more difficult for X-nuclei dynamic MRSI (without hyperpolarization) because of their relatively low sensitivities. This work presents a new method for high-resolution dynamic 31P-MRSI using high-order partially separable (HOPS) functions.²

We use a special time interleaved data acquisition scheme for sparse sampling (an example is shown in Fig.1). The scheme collects three complementary datasets: (1) a "training" dataset ($$$D_{1,f}$$$) acquired either before or after the dynamic process of interest, covering limited k-space but with a high spectral bandwidth, to determine the spectral distribution of the image function; (2) another "training" dataset ($$$D_{1,T}$$$) acquired during the dynamic process, covering limited k-space but at a high temporal sampling rate, to capture the dynamics of the image function; and (3) a sparsely sampled "imaging" dataset ($$$D_2$$$) also acquired during the dynamic process, interleaved with $$$D_{1,T}$$$, and covering extended k-space to recover the image function at a high resolution.

The image function $$$\rho(\mathbf{x},f,T)$$$ (representing spatial-spectral-temporal distribution) is reconstructed from the sparse data by taking advantage of data correlation in multiple dimensions. We represent $$$\rho(\mathbf{x},f,T)$$$ using HOPS functions: $$\rho(\mathbf{x},f,T)=\sum_{l=1}^{L} \sum_{m=1}^{M} \sum_{n=1}^{N} c_{lmn} \theta_{l}(\mathbf{x}) \phi_{m}(f) \psi_{n}(T)$$ where $$$\{\theta_{l}(\mathbf{x})\}_{l=1}^{L}$$$, $$$\{\phi_{m}(f)\}_{m=1}^{M}$$$, and $$$\{\psi_{n}(T)\}_{n=1}^{N}$$$ are basis functions that describe the variation of the image function along the spatial, spectral, and temporal axes, respectively, and $$$\{ c_{lmn} \}_{l,m,n=1}^{L,M,N}$$$ are the corresponding coefficients. Mathematically, this model leads to low-rank tensors, e.g., $$$\rho(\mathbf{x},f,T)$$$ can be represented as an order-3 tensor in the Tucker form³ with a tensor rank ($$$L$$$,$$$M$$$,$$$N$$$) after discretization.

We use an explicit subspace pursuit approach to reconstruction. First, the spectral and temporal basis functions ($$$\{\phi_{m}(f)\}_{m=1}^{M}$$$ and $$$\{\psi_{n}(T)\}_{n=1}^{N}$$$) are determined from the "training" dataset $$$D_{1,f}$$$ and $$$D_{1,T}$$$, respectively, using SVD.^2,4 The $$$B_0$$$ inhomogeneity effects on $$$D_{1,f}$$$ are corrected using the method in Ref. 5 prior to the estimation. Second, denoting the estimated basis functions as $$$\{\hat{\phi}_{m}(f)\}_{m=1}^{M}$$$ and $$$\{\hat{\psi}_{n}(T)\}_{n=1}^{N}$$$, $$$\rho(\mathbf{x},f,T)$$$ is recovered by determining $$$\{ c_{lmn} \}_{l,m,n=1}^{L,M,N}$$$ and $$$\{\theta_{l}(\mathbf{x})\}_{l=1}^{L}$$$ via fitting the data in $$$D_2$$$:$$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}(f) \hat{\psi}_{n}(T)\}\parallel_2^2 + R_{1}(\theta_{l}(\mathbf{x})) + R_{2}(c_{lmn}),$$ where $$$\mathbf{s}_{2}$$$ contains the $$$D_2$$$ data, $$$\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 using in vivo dynamic 31P-MRSI experiments (approved by our local IRB), which were carried out on a 3T Siemens Trio scanner equipped with a dual-channel 31P surface coil (PulseTeq, UK).

The data acquisition protocol is as follows. A low-resolution CSI acquisition for $$$D_{1,T}$$$ (TR/TE=160/3 ms, 2 kHz BW, 4x4x4 spatial encodings, and 2 averages) was interleaved with a high-resolution EPSI acquisition for $$$D_{2}$$$ (the same TR/TE, 40 kHz BW, 32x32x12 spatial encodings, 64 echoes, and bipolar acquisition). The interleaved acquisition was repeated 12 times, followed by another CSI acquisition for $$$D_{1,f}$$$ (the same TR/TE and BW, 8x8x6 spatial encodings, and 6 averages). We performed the entire acquisition two times. First the subject was asked to keep still while we collected a static dataset. The subject was then asked to perform repeated plantar flexion and dorsiflexion exercises without resistance for 5 minutes to stress the muscles of the lower leg. After exercise, the subject was again asked to remain still while the same protocol was used to collect a dynamic dataset to observe recovery.

Figure 2 shows a set of representative dynamic 31P-MRSI results obtained by the proposed method. The reconstructed PCr map and spectra had both high resolution (7.5 mm x 7.5 mm x 10 mm) and high SNR, as shown in Figs. 2b and 2c. Figure 2 plots the changes of the PCr peak over time, which were flat before exercise and showed an expected exponential recovery after the exercise.

This work presents a new approach to high-resolution dynamic MRSI using high-order partially separable functions. The proposed method has been validated using in vivo dynamic 31P-MRSI experiments, producing very encouraging results. It could enable a range of new applications of dynamic MRSI.