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.