High-Resolution Dynamic 31P-MRSI Using a Low-Rank Tensor Model

Chao Ma; Bryan Clifford; Yuchi Liu; Yuning Gu; Fan Lam; Xin Yu; Zhi-Pei Liang

doi:10.1002/mrm.26762

. Author manuscript; available in PMC: 2018 Aug 1.

Published in final edited form as: Magn Reson Med. 2017 May 28;78(2):419–428. doi: 10.1002/mrm.26762

High-Resolution Dynamic ³¹P-MRSI Using a Low-Rank Tensor Model

Chao Ma ^1,², Bryan Clifford ^2,³, Yuchi Liu ⁴, Yuning Gu ⁴, Fan Lam ^2,³, Xin Yu ⁴, Zhi-Pei Liang ^2,^3,^*

PMCID: PMC5562044 NIHMSID: NIHMS889879 PMID: 28556373

Abstract

Purpose

To develop a rapid ³¹P-MRSI method with high spatiospectral resolution using low-rank tensor-based data acquisition and image reconstruction.

Methods

The multidimensional image function of ³¹P-MRSI is represented by a low-rank tensor to capture the spatial–spectral–temporal correlations of data. A hybrid data acquisition scheme is used for sparse sampling, which consists of a set of “training” data with limited k-space coverage to capture the subspace structure of the image function, and a set of sparsely sampled “imaging” data for high-resolution image reconstruction. An explicit subspace pursuit approach is used for image reconstruction, which estimates the bases of the subspace from the “training” data and then reconstructs a high-resolution image function from the “imaging” data.

Results

We have validated the feasibility of the proposed method using phantom and in vivo studies on a 3T whole-body scanner and a 9.4T preclinical scanner. The proposed method produced high-resolution static ³¹P-MRSI images (i.e., 6.9×6.9×10 mm³ nominal resolution in a 15-min acquisition at 3T) and high-resolution, high-frame-rate dynamic ³¹P-MRSI images (i.e., 1.5 × 1.5 × 1.6 mm³ nominal resolution, 30 s/frame at 9.4T).

Conclusions

Keywords: 31P-MRSI, dynamic 31P-MRSI, partial separability, subspace modeling, low-rank matrix, low-rank tensor

INTRODUCTION

Phosphorus-31 magnetic resonance spectroscopy and imaging (³¹P-MRS/I) provides a unique capability for noninvasive quantification of the high-energy phosphate metabolites, which is useful for direct assessment of mitochondrial energy metabolism in vivo (1). In particular, monitoring the depletion and resynthesis of phosphocreatine (PCr) during exercise-recovery or ischemia-reperfusion by dynamic ³¹P-MRS/I allows for the assessment of mitochondrial oxidative capacity in the muscle (2,3). Furthermore, ³¹P magnetization-transfer provides quantitative evaluation of the synthesis of high-energy phosphate metabolites via creatine kinase and ATP synthase (4). However, because the concentrations of phosphate metabolites are at least three-orders of magnitude lower than water protons, current ³¹P-MRS methods require very long acquisition time to achieve adequate signal-to-noise ratio (SNR), significantly limiting its clinical utility. Despite the significant progress on fast imaging sequences (5–10) and advanced image reconstruction methods (11,12), the applications of ³¹P-MRSI are still limited by long data acquisition time, poor spatial resolution, and low SNR. To alleviate the limitations, most ³¹P-MRS studies have employed either non-localized or single voxel techniques, rendering the assessment of metabolic heterogeneity impossible.

Recently, a subspace based approach, called SPICE (SPectroscopic Imaging by exploiting spatiospectral CorrElation), has been proposed for accelerated high-resolution MRSI (13). In SPICE, MRSI signals are represented by a low-rank model (14–17), which captures the spatiospectral correlation of the data. This model implies that high-dimensional MRSI data reside in a very low dimensional subspace, thus enabling recovery of high resolution MRSI images from undersampled k-space data via special data acquisition schemes (e.g., sparse sampling) and image reconstruction methods. SPICE has been successfully applied to proton MRSI (¹H-MRSI) of the brain, producing 3D MRSI images at 3 mm isotropic resolution from a less than 10-min acquisition (18–20). Two unique properties of ³¹P-NMR make application of SPICE to ³¹P-MRSI promising. First, the ³¹P spectra of most detectable metabolites have a much simpler structure than the ¹H spectra, containing only a small number of resonance peaks with a large chemical shift dispersion of approximately 30 ppm. Second, while removing the nuisance water and lipid signals in ¹H-MRSI is very challenging (21), water and lipids do not contribute to ³¹P-MRSI signals. These properties make it possible to model ³¹P-MRSI signals using the low-rank model with a lower rank than ¹H-MRSI and therefore a reduced number of degrees of freedom.

This work presents a generalization of the SPICE method for both high-resolution static and dynamic ³¹PMRSI. We extend the low-rank matrix model in SPICE to a low-rank tensor model (14,22). This model allows us to take advantage of the inherent correlations of high-dimensional data in multiple dimensions, i.e., the spatial–spectral–temporal correlations of the dynamic ³¹PMRSI signal, which further reduces the number of degrees of freedom of the signal. A hybrid data acquisition scheme is used for sparse sampling, which consists of a set of “training” data with limited k-space coverage to capture the subspace structure (spectral and temporal basis functions) of the image function, and a set of sparse “imaging” data for high-resolution image reconstruction. An explicit subspace pursuit approach is used for image reconstruction, which estimates the bases of the subspace from the “training” data and then reconstructs a high-resolution image function by fitting the subspace model to the “imaging” data using the estimated bases.

The feasibility of the proposed method was demonstrated using phantom and in vivo studies on a 3T whole-body scanner and a 9.4T preclinical scanner. A preliminary account of this work was presented in (23–25).

THEORY

In the following, we describe the signal model, data acquisition scheme and reconstruction procedure of the proposed method in the context of dynamic ³¹P-MRSI and treat static ³¹P-MRSI as a special case.

Signal Model

Denote the image function of dynamic MRSI as ρ(x, f, T), where x, f, and T denote the spatial, spectral, and temporal axis, respectively. In practice, the variation of ρ(x, f, T) along each of these axes can be approximated by a linear combination of a small number of basis functions. For instance, the spectral distribution of ρ(x, f, T) can be represented by a finite number of resonance peaks (e.g., PCr, Pi, and the α, β, and γ peaks of ATP). The temporal variation of ρ(x, f, T) can be represented by a set of exponential functions with different decay or recovery time constants. In order to take advantage of these properties, we propose to represent the image function ρ(x, f, T) using a unified and more general model, known as the partially separable (PS) model (14):

ρ (x, f, T) = \sum_{l = 1}^{L} \sum_{m = 1}^{M} \sum_{n = 1}^{N} c_{l, m, n} θ_{l} (x) ϕ_{m} (f) ψ_{n} (T),

[1]

where ${θ_{l} (x)}_{l = 1}^{L}$ , ${ϕ_{m} (f)}_{m = 1}^{M}$ , and ${ψ_{n} (T)}_{n = 1}^{N}$ denote basis functions that describe the variation of ρ(x, f, T) along the spatial, spectral, and temporal axis, respectively, and ${c_{l, m, n}}_{l, m, n = 1}^{L, M, N}$ the corresponding model coefficients.

In the case of static MRSI, the model in Equation [1] is reduced to

ρ (x, f) = \sum_{l = 1}^{L} \sum_{m = 1}^{M} c_{l, m} θ_{l} (x) ϕ_{m} (f) .

[2]

Without loss of generality, assuming L ≤ M, Equation [2] can be rewritten as

ρ (x, f) = \sum_{l = 1}^{L} θ_{l} (x) (\sum_{m = 1}^{M} c_{l, m} ϕ_{m} (f)), = \sum_{l = 1}^{L} θ_{l} (x) {\tilde{ϕ}}_{l} (f),

[3]

where the model orders of the spatial and spectral basis functions become the same. The model in Equation [3] is the same as in our previous low-rank based methods for static MRSI (16,18–20).

Mathematically, the PS model in Equation [1] implies a low-rank tensor structure. More specifically, after discretization, the image function ρ(x, f, T) can be represented by a three-way array be or a tensor $ρ \in ℂ^{P \times Q \times R}$ , where P, Q, and R denote the dimension of ρ along the spatial, spectral, and temporal axis, respectively. For a given ρ(x, f, T), P is determined by the FOV and desired spatial resolution, Q by the spectral bandwidth and desired spectral resolution, and R by the length of the dynamic process of interest and desired frame rate. According to Equation [1], ρ can be expressed in the Tucker form (26) as follows:

ρ = \sum_{l = 1}^{L} \sum_{m = 1}^{M} \sum_{n = 1}^{N} c_{l, m, n} θ_{l} ○ ϕ_{m} ○ ψ_{n},

[4]

where θ_l, ϕ_m, and ψ_n are vectors that concatenate samples of θ_l(x), ϕ_m(f), and ψ_n(T) in Equation [1] on a discrete grid defined by ρ and ‘^○’ denotes the vector outer product (27). The rank of the tensor in Equation [4] is specified by (L, M, N), known as the multilinear rank of a tensor (27). In practice, L, M, and N are much smaller than P, Q, and R, respectively. This implies that the high-dimensional data ρ has a lot of redundancy (i.e., is highly correlated). Compared to the low-rank matrix model that exploits data correlation in two dimensions (e.g., the spatiospectral correlations in static MRSI data), the low-rank tensor model exploits data correlation in multiple dimensions simultaneously (e.g., the spatial–spectral–temporal correlations in dynamic MRSI data).

Data Acquisition

We propose a hybrid data acquisition scheme to exploit the PS property of the image function ρ(x, f, T) for high-resolution dynamic MRSI. The proposed scheme collects two complementary datasets: a “training” dataset for estimation of the subspace structure of ρ(x, f, T), and a sparsely sampled “imaging” dataset for high-resolution reconstruction of ρ(x, f, T). The “training” dataset further consists of two subsets, one of which (denoted as $D_{1, f}$ ) is acquired to estimate the spectral basis ( ${ϕ_{m}}_{m = 1}^{M}$ in Eq. [4]), and the other (denoted as $D_{1, T}$ ) to estimate the temporal basis ( ${ψ_{n}}_{n = 1}^{N}$ in Eq. [4]). $D_{1, f}$ can be acquired either before or after the dynamic process being studied with a high spectral bandwidth, using a chemical shift imaging (CSI) sequence or a fast echo-planar spectroscopic imaging (EPSI) sequence (if SNR allows). $D_{1, T}$ , on the other hand, is acquired at a high temporal rate during the dynamic process being studied along with the “imaging” dataset (denoted as $D_{2}$ ) either in an inter-leaved fashion or as part of $D_{2}$ . $D_{1, f}$ and $D_{1, T}$ only need to cover limited k-space for accurate subspace estimation. The “imaging” dataset $D_{2}$ covers an extended k-space to determine the spatial basis ( ${θ_{l}}_{l = 1}^{L}$ in Eq. [4]) at high resolution. $D_{2}$ can be sampled sparsely in the (k, t, T)-space. Many fast sequences such as EPSI sequences with Cartesian or spiral trajectories (6,7) can be used to acquire $D_{2}$ . An example of the proposed data acquisition scheme is illustrated in Figure 1.

FIG. 1 — An example of the proposed acquisition scheme. The “training” dataset $D_{1, T}$ is acquired at high temporal sampling rate but with limited k-space coverage (blue boxes). The sparse “imaging” dataset $D_{2}$ is acquired with extended k-space coverage (green cross signs). The acquisitions of $D_{1, T}$ and $D_{2}$ are interleaved in the whole dynamic process of the imaging object. Another “training” dataset $D_{1, f}$ is acquired with high spectral bandwidth but limited k-space coverage (red box) when the studied dynamic process is relatively stationary. Signal averaging can be used to improve SNR when acquiring $D_{1, f}$ (as reflected by the width of the red box).

Image Reconstruction

We adopt an explicit subspace pursuit approach for image reconstruction, which first estimates the bases of the subspace from the “training” data and then reconstructs the high-dimensional image by fitting the subspace model to the “imaging” data with the estimated bases.

Subspace Estimation

To estimate the spectral basis from $D_{1, f}$ , the B₀ inhomogeneity effects are first compensated using the method in (28), which aims to correct the B₀ inhomogeneity effects on MRSI data with limited k-space coverage by using a measured high-resolution B₀ inhomogeneity map and anatomical constraints, i.e., edge information derived from ¹H structural images. We then perform low-rank approximation (LORA) denoising to improve the SNR of $D_{1, f}$ (16), which takes advantage of both the PS and linear predictability properties of MRSI data for denoising. Denote the resultant data as ${\hat{D}}_{1, f} = {{\hat{d}}_{1, f} (k_{p}, f_{q})}_{p, q = 1}^{P_{1, f}, Q_{1, f}}$ , where P₁,_f and Q₁,_f are the number of samples along k and f axis, respectively. The spectral basis ${ϕ_{m}}_{m = 1}^{M}$ are estimated by calculating the M principal left-singular vectors of the following Casorati matrix:

C ({\hat{D}}_{1, f}) = [\begin{matrix} {\hat{d}}_{1, f} (k_{1}, f_{1}) & \dots & {\hat{d}}_{1, f} (k_{p, f_{1}}) & \dots & {\hat{d}}_{1, f} (k_{P_{1, f}}, f_{1}) \\ {\hat{d}}_{1, f} (k_{1}, f_{2}) & \dots & {\hat{d}}_{1, f} (k_{p}, f_{2}) & \dots & {\hat{d}}_{1, f} (k_{P_{1, f}}, f_{2}) \\ ⋮ & ⋱ & ⋮ & ⋱ & ⋮ \\ {\hat{d}}_{1, f} (k_{1}, f_{Q_{1, f}}) & \dots & {\hat{d}}_{1, f} (k_{p}, f_{Q_{1, f}}) & \dots & {\hat{d}}_{1, f} (k_{P_{1, f}}, f_{Q_{1, f}}) \end{matrix}]

[5]

See Refs. 13,28,29 for more details.

Assuming B₀ inhomogeneity does not change over time, it is straightforward to estimate the temporal basis from $D_{1, T}$ . Denote the k-space samples in $D_{1, T}$ as ${d_{1, T} (k_{p}, f_{q}, T_{r})}_{p, q, r = 1}^{P_{1, T}, Q_{1, T}, R}$ , where P₁,_T, Q₁,_T, and R are the number of samples along k, f, and T axis, respectively. We first perform LORA denoising on $D_{1, T}$ at each T_r. Denote the resultant data as ${\hat{D}}_{1, T}$ . Note that ${\hat{D}}_{1, T}$ itself is a tensor. In order to estimate the temporal basis functions, we unfold ${\hat{D}}_{1, T}$ and form the following Casorati matrix:

C ({\hat{D}}_{1, T}) = [\begin{matrix} {\hat{d}}_{1, T} (k_{1}, f_{1}, T_{1}) & \dots & {\hat{d}}_{1, T} (k_{P_{1}, T}, f_{1}, T_{1}) & \dots & {\hat{d}}_{1, T} (k_{p}, f_{q}, T_{1}) & \dots & {\hat{d}}_{1, T} (k_{P_{1}, T}, f_{Q_{1}, T}, T_{1}) \\ {\hat{d}}_{1, T} (k_{1}, f_{1}, T_{2}) & \dots & {\hat{d}}_{1, T} (k_{P_{1}, T}, f_{1}, T_{2}) & \dots & {\hat{d}}_{1, T} (k_{p}, f_{q}, T_{2}) & \dots & {\hat{d}}_{1, T} (k_{P_{1}, T}, f_{Q_{1}, T}, T_{2}) \\ ⋮ & ⋱ & ⋮ & ⋱ & ⋮ \\ {\hat{d}}_{1, T} (k_{1}, f_{1}, T_{R}) & \dots & {\hat{d}}_{1, T} (k_{P_{1}, T}, f_{1}, T_{R}) & \dots & {\hat{d}}_{1, T} (k_{p}, f_{q}, T_{R}) & \dots & {\hat{d}}_{1, T} (k_{P_{1}, T}, f_{Q_{1}, T}, T_{R}) \end{matrix}] .

[6]

The temporal basis ${ψ_{n}}_{n = 1}^{N}$ are then estimated by calculating the N principal left-singular vectors of the Casorati matrix in Equation [6].

Reconstruction

Denoting the estimated spectral and temporal bases as ${{\hat{ϕ}}_{m}}_{m = 1}^{M}$ and ${{\hat{ψ}}_{n}}_{n = 1}^{N}$ , ρ in Equation [4] is reconstructed by determining the model coefficients ${c_{l, m, n}}_{l, m, n = 1}^{L, M, N}$ and spatial basis ${θ_{l}}_{l = 1}^{L}$ via fitting the “imaging” data in $D_{2}$ (22):

min_{{c_{l, m, n}}_{l, m, n = 1}^{L, M, N}, {θ_{l}}_{l = 1}^{L}} {‖ d_{2} - Ω ℱ_{B} {\sum_{l = 1}^{L} \sum_{m = 1}^{M} \sum_{n = 1}^{N} c_{l, m, n} θ_{l} ○ {\hat{ϕ}}_{m} ○ {\hat{ψ}}_{n}} ‖}_{2}^{2} + R ({θ_{l}}_{l = 1}^{L}),

[7]

where d₂ denotes a vector concatenating k-space samples in $D_{2}$ , Ω a sampling pattern in the (k, f, T)-space, and $ℱ_{B}$ a Fourier encoding operator taking B₀ field inhomogeneity into account. In Equation [7], the first term is a standard data consistency penalty and the second term is used to incorporate the prior knowledge of the spatial distributions of metabolites (30).

While the optimization problem in Equation [7] can be solved by updating ${c_{l, m, n}}_{l, m, n = 1}^{L, M, N}$ and ${θ_{l}}_{l = 1}^{L}$ in an alternating scheme as described in Ref. 22, we propose to solve the problem using the following algorithm because of its computational efficiency:

Estimate an MRSI image (denoted as $ρ_{r}^{(1)} = {ρ_{p, q, r}^{(1)})}_{p, q = 1}^{P, Q}$ ) at each time point T_r, r = 1, …, R by solving the following optimization problem:
${{\hat{a}}_{m, r}}_{m = 1}^{M} = \underset{{a_{m, r}}_{m = 1}^{M}}{arg min} {‖ d_{2, r} - Ω_{r} ℱ_{B} {\sum_{m = 1}^{M} a_{m, r} ○ {\hat{ϕ}}_{m}} ‖}_{2}^{2} + R ({a_{m, r}}_{m = 1}^{M}),$ [8]
where d₂,_r denotes a vector concatenating the k-space samples in $D_{2}$ at T_r, Ω_r the sampling pattern in the (k, f, T)-space at T_r, and a_m,_r the corresponding spatial coefficients at T_r to be estimated.
Project each ${\hat{a}}_{m, r}$ r = 1, …, R onto a subspace spanned by the temporal subspace ${{\hat{ψ}}_{n}}_{n = 1}^{N}$ :
${\hat{b}}_{m, n} = \underset{{b_{m, n}}_{n = 1}^{N}}{arg min} {‖ C ({{\hat{a}}_{m, r}}_{r = 1}^{R}) - \sum_{n = 1}^{N} b_{m, n} ○ {\hat{ψ}}_{n} ‖}_{F}^{2}$ [9]
where $C ({{\hat{a}}_{m, r}}_{r = 1}^{R}) \in ℂ^{P \times R}$ is a Casorati matrix formed by ${\hat{a}}_{m, r}$ . Denote the resultant image as ρ⁽²⁾:
$ρ^{(2)} = \sum_{n = 1}^{N} \sum_{m = 1}^{M} {\hat{b}}_{m, n} ○ {\hat{ϕ}}_{m} ○ {\hat{ψ}}_{n} .$ [10]
Estimate the spatial basis functions by calculating the L principal left-singular vectors of the Casorati matrix formed by ρ⁽²⁾, i.e., $C (ρ^{(2)}) \in ℂ^{P \times Q R}$ . Denote the estimated spatial basis functions as ${{\hat{θ}}_{l}}_{l = 1}^{L}$ .
Estimate the model coefficients by solving the following least-squares problem:
${{\hat{c}}_{l, m, n}}_{l, m, n = 1}^{L, M, N} = \underset{{c_{l, m, n}}_{l, m, n = 1}^{L, M, N}}{arg min} {‖ ρ^{(2)} - \sum_{l = 1}^{L} \sum_{m = 1}^{M} \sum_{n = 1}^{N} c_{l, m, n} {\hat{θ}}_{l} ○ {\hat{ϕ}}_{m} ○ {\hat{ψ}}_{n} ‖}_{2}^{2} .$ [11]

METHODS

Static and Dynamic ³¹P-MRSI at 3T

The feasibility of the proposed method was demonstrated using phantom and in vivo studies (approved by our local IRB) on a 3T Siemens Trio scanner (Siemens Healthcare, Germany) equipped with a dual-channel ³¹P surface coil (PulseTeq, UK).

A static ³¹P-MRSI study was performed on a phantom to compare the performance of the proposed method with the conventional CSI and EPSI methods. The phantom consisted of one shorter and wider tube with an inner diameter of 25 mm and a length of 95 mm and two longer and thinner tubes with an inner diameter of 15 mm and a length of 125 mm. The three vials were filled with 100 mM sodium phosphate (Sigma Aldrich) solution and mounted in a cylindrical jar filled with NaCl-doped water. Experiments were carried out with acquisitions equivalent in time. 3D ³¹P-MRSI data were acquired using customized free induction decay (FID) MRSI sequences with 170 ms TR, 3.3 ms TE, 18° flip angle, and 240 × 240 × 120 mm³ field-of-view (FOV). The proposed method consisted of two acquisitions. The “training” dataset $D_{1, f}$ was acquired using a CSI sequence with 8 × 8 × 8 spatial encodings, 192 spectral encodings, 2 kHz sampling bandwidth, and four signal averages (5.8 min acquisition). The “imaging” dataset $D_{2}$ was acquired using an EPSI sequence with 50 × 50 × 20 spatial encodings, 128 echoes (bipolar acquisition), 68 kHz sampling bandwidth, 0.96 ms echo spacing (the distance between two neighboring echoes), and four signal averages (11.2 min acquisition). For comparison, conventional CSI and EPSI data were also acquired using the same CSI and EPSI sequence, respectively. The CSI data were acquired with 10 × 10 × 10 spatial encodings and six signal averages. The EPSI data were acquired with 50 × 50 × 20 spatial encodings and six signal averages. The acquisition time of the CSI and EPSI experiment was 17 min as well. Additional ¹H structural images and B₀ inhomogeneity maps were acquired using a dual-echo gradient recalled echo (GRE) sequence and the body coil of the scanner with an in-plane resolution of 1.9 mm and a slice thickness of 5 mm.

An in vivo study was performed on a healthy subject to demonstrate the feasibility of the proposed method for static ³¹P-MRSI. The proposed method was used to acquire 3D ³¹P-MRSI data from the calf muscle of the subject using a FID-CSI sequence. The “training” dataset $D_{1, f}$ was acquired with 170 ms TR, 3.0 ms TE, 18° flip angle, 220 × 220 × 120 mm³ FOV, 8 × 8 × 6 spatial encodings, 192 spectral encodings, 2 kHz sampling bandwidth, and two signal averages. The “imaging” dataset $D_{2}$ was acquired with the same TR, TE, FOV, and sampling bandwidth, 32 × 32 × 12 spatial encodings (6.9 × 6.9 × 10 mm³ nominal resolution), 192 spectral encodings, elliptical sampling, and single signal averages. The total acquisition of the proposed method was 15 min. Additional ¹H structural images and B₀ inhomogeneity maps were acquired using a dual-echo GRE sequence and the body coil of the scanner with an in-plane resolution of 1.9 mm and a slice thickness of 5 mm.

Another in vivo study was performed on a healthy subject to demonstrate the feasibility of the proposed method for dynamic ³¹P-MRSI. First, a set of $D_{1, f}$ , $D_{1, T}$ , and $D_{2}$ data were acquired while the subject kept still. $D_{1, f}$ data were acquired using an FID-CSI sequence with 160 ms TR, 3 ms TE, 240 × 240 × 120 mm³ FOV, 8 × 8 × 8 spatial encodings, 192 spectral encodings, 2 kHz sampling bandwidth, and six averages. $D_{1, T}$ and $D_{2}$ data were then acquired using an FID-CSI and FID-EPSI sequence, respectively, in an interleaved fashion. More specifically, at each frame, $D_{1, T}$ data were acquired using the FID-CSI sequence with the same TR, TE, and FOV, 4 × 4 × 4 spatial encodings, and two averages (22 s acquisition), followed by an EPSI acquisition for $D_{2}$ with the same TR, TE, and FOV, 32 × 32 × 12 spatial encodings (7.5 × 7.5 × 10 mm³ nominal resolution), 128 echoes, 0.93 ms echo spacing, 40 kHz sampling bandwidth, and single average (65 s acquisition). The nominal frame rate was thus 87 s. The interleaved acquisition was repeated 12 times. The subject was then asked to perform repeated plantarflexion and dorsiflexion exercises without resistance for 5 min to stress the muscles of the lower leg. After exercise, the subject was again asked to remain still while the same imaging protocol was used to acquire another set of $D_{1, T}$ and $D_{2}$ data to observe the recovery of PCr.

Dynamic ³¹P-MRSI at 9.4T

To further demonstrate the feasibility of the proposed method for dynamic ³¹P-MRSI, in vivo data were acquired from the hindlimb of an anesthetized rat on a 9.4T Bruker scanner equipped with a custom-built ³¹P saddle coil and a Bruker ¹H volume coil. First, an FID-CSI sequence was used to acquire a $D_{1, f}$ dataset at baseline for 8 min with 160 ms TR, 0.69 ms TE, 17° flip angle, 24 × 24 × 20 mm³ FOV, 8 × 8 × 6 spatial encodings, 256 spectral encodings, 6 kHz sampling bandwidth, and eight signal averages. Ischemia was then induced by inflating a cuff placed around the animal’s thigh (see Ref. 30 for more details). An FID-spiral-EPSI sequence was used to acquire a $D_{2}$ dataset during 10 min ischemia and 5 min reperfusion with the same TR, TE, and FOV, 16 × 16 × 12 matrix size (stack of uniform density spirals), 20 echoes, 3.96 ms echo-spacing, 111.1 kHz sampling bandwidth, and 15 averages. The nominal frame rate of the spiral-EPSI acquisition was 30 s/frame. No $D_{1, T}$ data were acquired. Temporal basis functions were estimated from $D_{2}$ . The total experiment time was approximately 23 min.

RESULTS

The reconstruction dimensions and model orders of each experiment are summarized in Supporting Table S1.

Figure 2 shows the results from the equivalent-time acquisition experiment on a phantom. Figure 2a shows the ¹H structural images of the phantom. Figure 2b–d show the maps of Pi obtained by the conventional CSI, EPSI, and the proposed method, respectively. Figure 2e–g show representative spectra from Figure 2b–d. Compared to the conventional CSI and EPSI method, the proposed method was able to recover the spatiospectral distribution of the chemical compound with both high resolution and high SNR. Most notably, the three vials are clearly seen in the map of Pi obtained by the proposed method (Fig. 2d) and the corresponding spectrum (Fig. 2g) has an SNR comparable to the CSI result (Fig. 2e). The signal inhomogeneities were due to transmit and receive B₁ inhomogeneities of the dual-channel ³¹P surface coil. Note that images from three representative slices are shown in Figure 2a–d. Please see Supporting Figure S1 for images from all the slices.

FIG. 2 — Static ³¹P-MRSI on a phantom. a: ¹H structural images. b–d: Maps of Pi obtained by the conventional CSI (10 × 10 × 10 spatial encodings and six averages), EPSI (50 × 50 × 20 spatial encodings and six averages), and the proposed method (50 × 50 × 20 spatial encodings), respectively, in an equivalent-time acquisition experiment. **e–g**: Representative spectra from (b–d). The location of the spectrum is indicated by the red arrow in (a). Note that images from three representative slices are shown in (a–d). Please see Supporting Figure S1 for images from all the slices.

Figure 3 shows the static ³¹P-MRSI results of the calf muscle obtained by the proposed method on a healthy human subject. Figure 3a shows the ¹H structural images. Figure 3b,c show the maps of PCr obtained by taking spectral integral around the single peak of PCr. Figure 3d,e show representative spectra from Figure 3b,c, respectively. Compared to the results obtained by direct Fourier reconstruction of $D_{2}$ (Fig. 3b,d), the proposed method significantly improved the SNR (Fig. 3c,e). Besides the strong PCr peak, the peaks of Pi and ATP are clearly seen in the reconstructed spectrum and the bone areas with expected low ³¹P signals in the reconstructed PCr map that matched well with the ¹H structure image. Compared to the spectrum obtained by a low-resolution CSI acquisition with the same flip angle, TR and TE (Supporting Fig. S2), the Pi and ATP peaks in Figure 3e are at similar amplitudes. Also see Supporting Figure S3 for the five spectral basis functions that were used to produce the results in Figure 3.

FIG. 3 — Static ³¹P-MRSI of the calf muscle. a: ¹H structural images. b, c: Maps of PCr obtained by direct Fourier reconstruction of $D_{2}$ and the proposed method, respectively. d, e: Representative spectra from (b, c). The location of the spectrum is indicated by the red circle in (c).

Figures 4 and 5 together show the spatial–spectral–temporal distributions obtained by the proposed method in the dynamic ³¹P-MRSI experiment on a healthy subject. Figure 4 shows the results of the measurement before the subject performed repeated plantarflexion and dorsiflexion exercises. The time-average PCr map (Fig. 4b) and representative spectrum (Fig. 4c) show that the proposed method reconstructed the spatial–spectral distribution of the imaging object at both high spatial resolution and high SNR. Note that the results in Figure 3 were obtained using a CSI sequence, while those in Figure 4 were obtained using a hybrid CSI/EPSI sequence that has higher acquisition efficiency and is thus more suitable for a dynamic ³¹P-MRSI experiment. The representative changes of PCr over time at two voxels from different muscles groups are plotted as blue solid lines in Figure 5c,d, respectively, which show a constant PCr level as expected. After exercises, the same measurement was repeated. The representative changes of PCr over time at the same locations are plotted as red dashed lines in Figure 5c,d, respectively, showing expected exponential recovery after the exercise and distinct responses of different muscle groups.

FIG. 4 — Dynamic ³¹P-MRSI of the calf muscle. a: ¹H structural images. b: Time-average PCr map. c: Representative spectrum from (b). Its location is indicated by the red triangle in (b).

FIG. 5 — Dynamics of PCr from a representative slice. a: ¹H structural image. b: Time-average PCr map fused with the structural image in (a). c, d: Changes of PCr over time at two voxels from different muscle groups. The locations of the spectra in (c, d) are indicated by the black triangle and square in (b), respectively. The blue solid lines correspond to the measurement before the subject performed repeated plantarflexion and dorsiflexion exercises. The red dashed lines correspond to the measurement after the exercise.

Figures 6 and 7 show the dynamic ³¹P-MRSI results obtained by the proposed method in the animal experiment carried on a 9.4T preclinical system. Figure 6a shows the center five slices of the anatomical reference images of a rat leg. The corresponding PCr maps during ischemia (Fig. 6b–d) and reperfusion (Fig. 6e) show that PCr pool was significantly reduced during ischemia and rapidly replenished upon reperfusion. Dynamic spectra from two pixels in different muscle groups showed variation in the kinetics during ischemia-reperfusion (Fig. 7). These preliminary results demonstrate that our method was able to capture both the temporal dynamics and spatial variation of the image while maintaining good spectral quality. The dynamic behavior of the PCr spatiospectral distribution is consistent with the previous experiments reported in the literature (31–33). Note that the presented animal experiment was designed for proof-of-concept validation of our method. The SNR of the experiment can be further improved by using surface coils. The spatial resolution of the experiment can be improved by using a multishot variable density spiral trajectory to acquire $D_{1, T}$ and $D_{2}$ data simultaneously.

FIG. 6 — Slices of the anatomical reference images (a) and PCr peak integral (b–e) of the reconstruction at various time points. The pressure cuff was inflated and deflated at approximately 0 and 10 min respectively.

FIG. 7 — Spectra from two different spatial points in the reconstruction shown over time. The locations of the spectra in (a, b) correspond to the locations of the red circle and blue square in Figure 6b respectively.

DISCUSSION

We have proposed a low-rank tensor-based method for high-resolution ³¹P-MRSI, which is characterized by a hybrid data acquisition scheme and an explicit subspace pursuit approach to image reconstruction. According to the literature in matrix completion and recovery, it is possible to accurately estimate the subspace (row space) of a low-rank matrix (e.g., a rank-L matrix $A \in ℂ^{P \times Q}$ ) from a small number (e.g., s > L but s ≪ P) of fully sampled rows (34). The same applies to a low-rank tensor because a low-rank tensor can always be unfolded into a low-rank matrix. Therefore, the “training” datasets $D_{1, f}$ and $D_{1, T}$ only need to contain a small number of fully sampled data points in the (k,f)-plane and (k,T)-plane, respectively, for estimation of the tensor subspace. We choose to let $D_{1, f}$ and $D_{1, T}$ cover the center of k-space for SNR benefits. It can be shown that any spectral or temporal basis functions that are not captured by such “training” datasets are associated with spectral or temporal components that only appear as high-frequency spatial variations, which are rare for most applications. In practice, more samples are needed for $D_{1, f}$ to ensure accurate estimation of the spectral basis in the presence of B₀ field inhomogeneity. Based on our experience, 8 × 8 × 6 spatial encodings in $D_{1, f}$ are generally sufficient for typical 3D ³¹P-MRSI at 3T, which is still relatively small compared to the dimension of the target high-resolution reconstruction (e.g., 32×32×12 spatial encodings).

In image reconstruction, rank selection is critical. In principle, the rank of the proposed model is determined by the number of distinct metabolites, tissue types, and dynamic processes of the imaging object. In practice, we determine the rank adaptively by performing singular value decomposition (SVD) on the Casorati matrices formed by $D_{1, f}$ , $D_{1, T}$ , and ρ⁽²⁾ in Equation [10] and thresholding the calculated singular values. Automatic selection of rank is under investigation. Once the temporal and spectral basis functions are determined, the image reconstruction problem is reduced to estimating the model coefficients ${c_{l, m, n}}_{l, m, n = 1}^{L, M, N}$ and spatial basis ${θ_{l}}_{l = 1}^{L}$ via fitting the “imaging” data in $D_{2}$ as in Equation [7]. We proposed to estimate ${θ_{l}}_{l = 1}^{L}$ and ${c_{l, m, n}}_{l, m, n = 1}^{L, M, N}$ separately, i.e., via Equation [8] to [11], which works well when there are sufficient samples in the (k,f)-plane at each time frame of $D_{2}$ to allow successful frame-by-frame low-rank based reconstruction via Equation [8]. Our feasibility studies fell into this scenario. If a more general (k,f,T)-space sampling pattern is used to acquire $D_{2}$ , e.g., for higher frame rate, the joint estimation of ${θ_{l}}_{l = 1}^{L}$ and ${c_{l, m, n}}_{l, m, n = 1}^{L, M, N}$ in (22) can be used for image reconstruction, where the proposed method can be used to obtain the initial guess of ${θ_{l}}_{l = 1}^{L}$ and ${c_{l, m, n}}_{l, m, n = 1}^{L, M, N}$ .

Compared to the conventional CSI and EPSI method, our method has significant advantage in SNR as shown in Figures 2 and 3. Conceptually, the SNR benefit of our method comes from two sources. First, the spectral and temporal basis functions are estimated by performing SVD on the “training” data, which cover the center of the k-space and are often acquired with good SNR. Second, once the spectral and temporal basis functions are determined, the image reconstruction problem is reduced to determining the spatial basis functions and model coefficients from the “imaging” data. This significantly reduces the degrees-of-freedom of the imaging function, making it possible to recover the imaging function from the noisy “imaging” data with good SNR. Assuming no modeling error, the SNR of the low-rank tensor based reconstruction can be theoretically characterized by analyzing the noise properties of the estimated spectral, temporal, and spatial basis functions, respectively. The noise level of the spectral and temporal basis functions are determined by the noise level of the “training” data. The noise level of the estimated spatial basis functions and model coefficients can be characterized by performing perturbation analysis to the solution of the optimization problem in Equation [7]. This is our ongoing research and will be reported in future publications.

The general framework of the low-rank tensor based approach to high-dimensional MRI has been proposed and validated using proof-of-concept simulation studies in (22). This work presents a specialized data acquisition scheme (as shown in Fig. 1) and image reconstruction method (Eqs. [8] to [11]) for use with the low-rank tensor model to achieve high spatial-resolution, high frame-rate ³¹P-MRSI. Furthermore, the feasibility of the method in this work has been validated using phantom and in vivo experimental studies on a 3T scanner and a 9.4T preclincal scanner. Other applications of the low-rank tensor model include diffusion tensor imaging (22), parameter mapping (22), J-resolved MRSI (22,35), cardiac imaging (36), and dynamic electronic paramagnetic resonance imaging (37).

CONCLUSIONS

Dynamic spatiospectral variations of ³¹P-MRSI signals can be efficiently represented by a low-rank tensor. Exploiting this mathematical structure for data acquisition and image reconstruction can lead to fast ³¹P-MRSI with high resolution, frame-rate, and SNR. The feasibility of the proposed method was demonstrated using phantom and in vivo experiments on a healthy subject and on a rat undergoing an ischemia-reperfusion procedure. The proposed method could enable a range of new applications of ³¹P-MRSI.

Supplementary Material

Supporting Info

Table S1. Reconstruction dimension and model order.

Fig. S1. Static 31P-MRSI on a phantom. (a) 1H structural images. (b) to (d) Maps of Pi obtained by the conventional CSI (10 × 10 × 10 spatial encodings and 6 averages), EPSI (50 × 50 × 20 spatial encodings and 6 averages), and the proposed method (50 × 50 × 20 spatial encodings), respectively, in an equivalent-time acquisition experiment. (e) to (g) Representative spectra from (b) to (d). The location of the spectrum is indicated by the red arrow in (a).

Fig. S2. Static 31P-MRSI of the calf muscle by a low-resolution CSI acquisition (10 × 10 × 6 spatial encodings). (a) Map of PCr. (b) and (c) Representative spectrum before and after LORA (Nguyen et al., IEEE Trans Biomed Eng, 2013;60:78-89) denoising.

Fig. S3. Singular value distribution (a) and the cooresponding spectral basis functions that were used to generate Fig.3.

NIHMS889879-supplement-Supporting_Info.pdf^{(497.1KB, pdf)}

Acknowledgments

The authors thanked Dr. Ryan Larsen at the Beckman Institute at the University of Illinois at Urbana-Champaign for his assistance in the setup of ³¹P surface coil and preparation of ³¹P phantom.

Grant sponsor: National Institutes of Health; Grant numbers: NIH-1RO1-EB013695, NIH-R21EB021013-01, and NIH-R21-HL126215; Grant sponsor: Beckman Postdoctoral Fellowship (to C. Ma and F. Lam); Grant number: NIH-R01-CA165221 and NIH-R21-EB021710.

Footnotes

SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this article.

References

1.Chance B, Eleff S, Leigh JS. Noninvasive, nondestructive approaches to cell bioenergetics. Proc Natl Acad Sci. 1980;77:7430–7434. doi: 10.1073/pnas.77.12.7430. [DOI] [PMC free article] [PubMed] [Google Scholar]
2.Meyer RA. A linear model of muscle respiration explains monoexponential phosphocreatine changes. Am J Physiol Cell Physiol. 1988;254:C548–C553. doi: 10.1152/ajpcell.1988.254.4.C548. [DOI] [PubMed] [Google Scholar]
3.Paganini AT, Foley JM, Meyer RA. Linear dependence of muscle phosphocreatine kinetics on oxidative capacity. Am J Physiol. 1997;272:C501–C510. doi: 10.1152/ajpcell.1997.272.2.C501. [DOI] [PubMed] [Google Scholar]
4.Befroy DE, Rothman DL, Petersen KF, Shulman GI. 31P-magnetization transfer magnetic resonance spectroscopy measurements of in vivo metabolism. Diabetes. 2012;61:2669–2678. doi: 10.2337/db12-0558. [DOI] [PMC free article] [PubMed] [Google Scholar]
5.Pohmann R, von Kienlin M, Haase A. Theoretical evaluation and comparison of fast chemical shift imaging methods. J Magn Reson. 1997;129:145–160. doi: 10.1006/jmre.1997.1245. [DOI] [PubMed] [Google Scholar]
6.Posse S, Tedeschi G, Risinger R, Ogg R, Bihan DL. High speed 1H spectroscopic imaging in human brain by echo planar spatial-spectral encoding. Magn Reson Med. 1995;33:34–40. doi: 10.1002/mrm.1910330106. [DOI] [PubMed] [Google Scholar]
7.Adalsteinsson E, Irarrazabal P, Topp S, Meyer C, Macovski A, Spielman DM. Volumetric spectroscopic imaging with spiral-based k-space trajectories. Magn Reson Med. 1998;39:889–898. doi: 10.1002/mrm.1910390606. [DOI] [PubMed] [Google Scholar]
8.van der Kemp WJM, Boer VO, Luijten PR, Stehouwer BL, Veldhuis WB, Klomp DWJ. Adiabatic multi-echo 31P spectroscopic imaging (AMESING) at 7 T for the measurement of transverse relaxation times and regaining of sensitivity in tissues with short T2* values. NMR Biomed. 2013;26:1299–1307. doi: 10.1002/nbm.2952. [DOI] [PubMed] [Google Scholar]
9.Diaz AS, Akbari A, Noseworthy M. Proceedings of the ISMRM. Singapore: 2016. 31P MRSI of human calf muscles using flyback echo planar spectroscopic imaging (EPSI) readout gradients; p. 3936. [Google Scholar]
10.Valkovic L, Chmelik M, Meyerspeer M, Gagoski B, Rodgers CT, Krssak M, Andronesi OC, Trattnig S, Bogner W. Dynamic 31PMRSI using spiral spectroscopic imaging can map mitochondrial capacity in muscles of the human calf during plantar flexion exercise at 7 T. NMR Biomed. 2016;29:1825–1834. doi: 10.1002/nbm.3662. [DOI] [PMC free article] [PubMed] [Google Scholar]
11.Raghavan RS, Panda A, Valette J, James J, Heberlein K, Boettcher U, Henry P, Bansal N, Dydak U. Proceedings of the ISMRM. Honolulu, USA: 2009. 31P spectroscopic imaging with GRAPPA; p. 4317. [Google Scholar]
12.Hatay GH, Gumus C, Hakyemez B, Ozturk-Isik E. Proceedings of the ISMRM. Salt Lake City, USA: 2013. Accelerated phosphorus MR spectroscopic imaging of human brain using compressed sensing; p. 2017. [Google Scholar]
13.Lam F, Liang ZP. A subspace approach to high-resolution spectroscopic imaging. Magn Reson Med. 2014;71:1349–1357. doi: 10.1002/mrm.25168. [DOI] [PMC free article] [PubMed] [Google Scholar]
14.Liang ZP. Proceedings of IEEE International Symposium on Biomedical Imaging. USA: 2007. Spatiotemporal imaging with partially separable functions; pp. 988–991. [Google Scholar]
15.Haldar JP, Liang ZP. Proceedings of IEEE International Symposium on Biomedical Imaging. Rotterdam, Netherlands: 2010. Spatiotemporal imaging with partially separable functions: A matrix recovery approach; pp. 716–719. [Google Scholar]
16.Nguyen HM, Peng X, Do MN, Liang ZP. Denoising MR spectroscopic imaging data with low-rank approximations. IEEE Trans Biomed Eng. 2013;60:78–89. doi: 10.1109/TBME.2012.2223466. [DOI] [PMC free article] [PubMed] [Google Scholar]
17.Kasten J, Lazeyras F, Van De Ville D. Data-driven MRSI spectral localization via low-rank component analysis. IEEE Trans Med Imaging. 2013;32:1853–1863. doi: 10.1109/TMI.2013.2266259. [DOI] [PubMed] [Google Scholar]
18.Lam F, Ma C, Clifford B, Johnson CL, Liang ZP. High-resolution 1H-MRSI of the brain using SPICE: Data acquisition and image reconstruction. Magn Reson Med. 2016;76:1059–1070. doi: 10.1002/mrm.26019. [DOI] [PMC free article] [PubMed] [Google Scholar]
19.Ma C, Lam F, Qiang N, Johnson CL, Liang ZP. High-resolution 1H-MRSI of the brain using short-TE SPICE. Magn Reson Med. 2017;77:467–479. doi: 10.1002/mrm.26130. [DOI] [PMC free article] [PubMed] [Google Scholar]
20.Lam F, Lu H, Yang Y, Clifford B, Ma C, Robinson GE, Liang ZP. Proceedings of the ISMRM. Singapore: 2016. Ultrahigh-resolution metabolic imaging at 9.4 tesla; p. 385. [Google Scholar]
21.Ma C, Lam F, Johnson CL, Liang ZP. Removal of nuisance signals from limited and sparse 1H MRSI data using a union-of-subspaces model. Magn Reson Med. 2016;75:488–497. doi: 10.1002/mrm.25635. [DOI] [PMC free article] [PubMed] [Google Scholar]
22.He J, Liu Q, Christodoulou AG, Ma C, Lam F, Liang ZP. Accelerated high-dimensional MR imaging with sparse sampling using low-rank tensors. IEEE Trans Med Imaging. 2016;35:2119–2129. doi: 10.1109/TMI.2016.2550204. [DOI] [PMC free article] [PubMed] [Google Scholar]
23.Ma C, Lam F, Ning Q, Clifford B, Larsen R, Liang ZP. Proceedings of the ISMRM. Singapore: 2016. A subspace-based approach to high-resolution 31P-MRSI; p. 3956. [Google Scholar]
24.Ma C, Lam F, Ning Q, Clifford B, Liu Q, Johnson CL, Liang ZP. Proceedings of the ISMRM. Singapore: 2016. High-resolution dynamic 31P-MRSI using high-order partially separable functions; p. 875. [Google Scholar]
25.Liu Y, Clifford B, Ma C, Lam F, Liang ZP, Yu X. Proceedings of ISMRM. Honolulu: 2017. High resolution dynamic 31P-MRSI of ischemia-reperfusion in rat hindlimb at 9.4T using SPICE; p. 2951. [Google Scholar]
26.Tucker LR. Some mathematical notes on three-mode factor analysis. Psychometrika. 1966;31:279–311. doi: 10.1007/BF02289464. [DOI] [PubMed] [Google Scholar]
27.Kolda TG, Bader BW. Tensor decompositions and applications. SIAM Rev. 2009;51:455–500. [Google Scholar]
28.Peng X, Nguyen H, Haldar JP, Hernando D, Wang XP, Liang ZP. Proceedings of the 32nd Annual Meeting of IEEE EMBC. Argentina: 2010. Correction of field inhomogeneity effects on limited k-space MRSI data using anatomical constraints; pp. 883–886. [DOI] [PubMed] [Google Scholar]
29.Liu Y, Ma C, Clifford B, Lam F, Curtis LJ, Liang ZP. Improved low-rank filtering of magnetic resonance spectroscopic imaging data corrupted by noise and B0 field inhomogeneity. IEEE Trans Biomed Eng. 2016;63:841–849. doi: 10.1109/TBME.2015.2476499. [DOI] [PubMed] [Google Scholar]
30.Haldar JP, Hernando D, Song SK, Liang ZP. Anatomically constrained reconstruction from noisy data. Magn Reson Med. 2008;59:810–818. doi: 10.1002/mrm.21536. [DOI] [PubMed] [Google Scholar]
31.Liu YC, Mei XB, Li JL, Lai N, Yu X. Mitochondrial function assessed by 31P MRS and BOLD MRI in non-obese type 2 diabetic rats. Physiol Rep. 2016;4:e12890. doi: 10.14814/phy2.12890. [DOI] [PMC free article] [PubMed] [Google Scholar]
32.Morikawa S, Kido C, Inubushi T. Observation of rat hind limb skeletal muscle during arterial occlusion and reperfusion by 31P MRS and 1H MRI. Magn Reason Imaging. 1991;9:269–274. doi: 10.1016/0730-725x(91)90411-e. [DOI] [PubMed] [Google Scholar]
33.Morikawa S, Inubushi T, Kito C. Heterogeneous metabolic changes in the calf muscle of the rat during ischemia-reperfusion: in vivo analysis by 31P nuclear magnetic resonance chemical shift imaging and 1H magnetic resonance imaging. Cardiovasc Surg. 1993;1:337–342. [PubMed] [Google Scholar]
34.Guruswami V, Sinop AK. Proceedings of the ACM-SIAM Symposium on Discrete Algorithms. Kyoto, Japan: 2012. Optimal column-based low-rank matrix reconstruction; pp. 1207–1214. [Google Scholar]
35.Ma C, Lam F, Liu Q, Liang ZP. Proceedings of the ISMRM. Singapore: 2016. Accelerated high-resolution multidimensional 1H-MRSI using low-rank tensors; p. 379. [Google Scholar]
36.Christodoulou AG, Shaw JL, Sharif B, Li D. Proceedings of the ISMRM. Singapore: 2016. A general low-rank tensor framework for high-dimensional cardiac imaging: Application to time-resolved T1 mapping; p. 867. [Google Scholar]
37.Christodoulou AG, Redler G, Clifford B, Liang ZP, Halpern HJ, Epel B. Fast dynamic electron paramagnetic resonance (EPR) oxygen imaging using low-rank tensors. J Magn Reson. 2016;270:176–182. doi: 10.1016/j.jmr.2016.07.006. [DOI] [PMC free article] [PubMed] [Google Scholar]

Associated Data

This section collects any data citations, data availability statements, or supplementary materials included in this article.

Supplementary Materials

Supporting Info

Table S1. Reconstruction dimension and model order.

Fig. S3. Singular value distribution (a) and the cooresponding spectral basis functions that were used to generate Fig.3.

NIHMS889879-supplement-Supporting_Info.pdf^{(497.1KB, pdf)}

[R1] 1.Chance B, Eleff S, Leigh JS. Noninvasive, nondestructive approaches to cell bioenergetics. Proc Natl Acad Sci. 1980;77:7430–7434. doi: 10.1073/pnas.77.12.7430. [DOI] [PMC free article] [PubMed] [Google Scholar]

[R2] 2.Meyer RA. A linear model of muscle respiration explains monoexponential phosphocreatine changes. Am J Physiol Cell Physiol. 1988;254:C548–C553. doi: 10.1152/ajpcell.1988.254.4.C548. [DOI] [PubMed] [Google Scholar]

[R3] 3.Paganini AT, Foley JM, Meyer RA. Linear dependence of muscle phosphocreatine kinetics on oxidative capacity. Am J Physiol. 1997;272:C501–C510. doi: 10.1152/ajpcell.1997.272.2.C501. [DOI] [PubMed] [Google Scholar]

[R4] 4.Befroy DE, Rothman DL, Petersen KF, Shulman GI. 31P-magnetization transfer magnetic resonance spectroscopy measurements of in vivo metabolism. Diabetes. 2012;61:2669–2678. doi: 10.2337/db12-0558. [DOI] [PMC free article] [PubMed] [Google Scholar]

[R5] 5.Pohmann R, von Kienlin M, Haase A. Theoretical evaluation and comparison of fast chemical shift imaging methods. J Magn Reson. 1997;129:145–160. doi: 10.1006/jmre.1997.1245. [DOI] [PubMed] [Google Scholar]

[R6] 6.Posse S, Tedeschi G, Risinger R, Ogg R, Bihan DL. High speed 1H spectroscopic imaging in human brain by echo planar spatial-spectral encoding. Magn Reson Med. 1995;33:34–40. doi: 10.1002/mrm.1910330106. [DOI] [PubMed] [Google Scholar]

[R7] 7.Adalsteinsson E, Irarrazabal P, Topp S, Meyer C, Macovski A, Spielman DM. Volumetric spectroscopic imaging with spiral-based k-space trajectories. Magn Reson Med. 1998;39:889–898. doi: 10.1002/mrm.1910390606. [DOI] [PubMed] [Google Scholar]

[R8] 8.van der Kemp WJM, Boer VO, Luijten PR, Stehouwer BL, Veldhuis WB, Klomp DWJ. Adiabatic multi-echo 31P spectroscopic imaging (AMESING) at 7 T for the measurement of transverse relaxation times and regaining of sensitivity in tissues with short T2* values. NMR Biomed. 2013;26:1299–1307. doi: 10.1002/nbm.2952. [DOI] [PubMed] [Google Scholar]

[R9] 9.Diaz AS, Akbari A, Noseworthy M. Proceedings of the ISMRM. Singapore: 2016. 31P MRSI of human calf muscles using flyback echo planar spectroscopic imaging (EPSI) readout gradients; p. 3936. [Google Scholar]

[R10] 10.Valkovic L, Chmelik M, Meyerspeer M, Gagoski B, Rodgers CT, Krssak M, Andronesi OC, Trattnig S, Bogner W. Dynamic 31PMRSI using spiral spectroscopic imaging can map mitochondrial capacity in muscles of the human calf during plantar flexion exercise at 7 T. NMR Biomed. 2016;29:1825–1834. doi: 10.1002/nbm.3662. [DOI] [PMC free article] [PubMed] [Google Scholar]

[R11] 11.Raghavan RS, Panda A, Valette J, James J, Heberlein K, Boettcher U, Henry P, Bansal N, Dydak U. Proceedings of the ISMRM. Honolulu, USA: 2009. 31P spectroscopic imaging with GRAPPA; p. 4317. [Google Scholar]

[R12] 12.Hatay GH, Gumus C, Hakyemez B, Ozturk-Isik E. Proceedings of the ISMRM. Salt Lake City, USA: 2013. Accelerated phosphorus MR spectroscopic imaging of human brain using compressed sensing; p. 2017. [Google Scholar]

[R13] 13.Lam F, Liang ZP. A subspace approach to high-resolution spectroscopic imaging. Magn Reson Med. 2014;71:1349–1357. doi: 10.1002/mrm.25168. [DOI] [PMC free article] [PubMed] [Google Scholar]

[R14] 14.Liang ZP. Proceedings of IEEE International Symposium on Biomedical Imaging. USA: 2007. Spatiotemporal imaging with partially separable functions; pp. 988–991. [Google Scholar]

[R15] 15.Haldar JP, Liang ZP. Proceedings of IEEE International Symposium on Biomedical Imaging. Rotterdam, Netherlands: 2010. Spatiotemporal imaging with partially separable functions: A matrix recovery approach; pp. 716–719. [Google Scholar]

[R16] 16.Nguyen HM, Peng X, Do MN, Liang ZP. Denoising MR spectroscopic imaging data with low-rank approximations. IEEE Trans Biomed Eng. 2013;60:78–89. doi: 10.1109/TBME.2012.2223466. [DOI] [PMC free article] [PubMed] [Google Scholar]

[R17] 17.Kasten J, Lazeyras F, Van De Ville D. Data-driven MRSI spectral localization via low-rank component analysis. IEEE Trans Med Imaging. 2013;32:1853–1863. doi: 10.1109/TMI.2013.2266259. [DOI] [PubMed] [Google Scholar]

[R18] 18.Lam F, Ma C, Clifford B, Johnson CL, Liang ZP. High-resolution 1H-MRSI of the brain using SPICE: Data acquisition and image reconstruction. Magn Reson Med. 2016;76:1059–1070. doi: 10.1002/mrm.26019. [DOI] [PMC free article] [PubMed] [Google Scholar]

[R19] 19.Ma C, Lam F, Qiang N, Johnson CL, Liang ZP. High-resolution 1H-MRSI of the brain using short-TE SPICE. Magn Reson Med. 2017;77:467–479. doi: 10.1002/mrm.26130. [DOI] [PMC free article] [PubMed] [Google Scholar]

[R20] 20.Lam F, Lu H, Yang Y, Clifford B, Ma C, Robinson GE, Liang ZP. Proceedings of the ISMRM. Singapore: 2016. Ultrahigh-resolution metabolic imaging at 9.4 tesla; p. 385. [Google Scholar]

[R21] 21.Ma C, Lam F, Johnson CL, Liang ZP. Removal of nuisance signals from limited and sparse 1H MRSI data using a union-of-subspaces model. Magn Reson Med. 2016;75:488–497. doi: 10.1002/mrm.25635. [DOI] [PMC free article] [PubMed] [Google Scholar]

[R22] 22.He J, Liu Q, Christodoulou AG, Ma C, Lam F, Liang ZP. Accelerated high-dimensional MR imaging with sparse sampling using low-rank tensors. IEEE Trans Med Imaging. 2016;35:2119–2129. doi: 10.1109/TMI.2016.2550204. [DOI] [PMC free article] [PubMed] [Google Scholar]

[R23] 23.Ma C, Lam F, Ning Q, Clifford B, Larsen R, Liang ZP. Proceedings of the ISMRM. Singapore: 2016. A subspace-based approach to high-resolution 31P-MRSI; p. 3956. [Google Scholar]

[R24] 24.Ma C, Lam F, Ning Q, Clifford B, Liu Q, Johnson CL, Liang ZP. Proceedings of the ISMRM. Singapore: 2016. High-resolution dynamic 31P-MRSI using high-order partially separable functions; p. 875. [Google Scholar]

[R25] 25.Liu Y, Clifford B, Ma C, Lam F, Liang ZP, Yu X. Proceedings of ISMRM. Honolulu: 2017. High resolution dynamic 31P-MRSI of ischemia-reperfusion in rat hindlimb at 9.4T using SPICE; p. 2951. [Google Scholar]

[R26] 26.Tucker LR. Some mathematical notes on three-mode factor analysis. Psychometrika. 1966;31:279–311. doi: 10.1007/BF02289464. [DOI] [PubMed] [Google Scholar]

[R27] 27.Kolda TG, Bader BW. Tensor decompositions and applications. SIAM Rev. 2009;51:455–500. [Google Scholar]

[R28] 28.Peng X, Nguyen H, Haldar JP, Hernando D, Wang XP, Liang ZP. Proceedings of the 32nd Annual Meeting of IEEE EMBC. Argentina: 2010. Correction of field inhomogeneity effects on limited k-space MRSI data using anatomical constraints; pp. 883–886. [DOI] [PubMed] [Google Scholar]

[R29] 29.Liu Y, Ma C, Clifford B, Lam F, Curtis LJ, Liang ZP. Improved low-rank filtering of magnetic resonance spectroscopic imaging data corrupted by noise and B0 field inhomogeneity. IEEE Trans Biomed Eng. 2016;63:841–849. doi: 10.1109/TBME.2015.2476499. [DOI] [PubMed] [Google Scholar]

[R30] 30.Haldar JP, Hernando D, Song SK, Liang ZP. Anatomically constrained reconstruction from noisy data. Magn Reson Med. 2008;59:810–818. doi: 10.1002/mrm.21536. [DOI] [PubMed] [Google Scholar]

[R31] 31.Liu YC, Mei XB, Li JL, Lai N, Yu X. Mitochondrial function assessed by 31P MRS and BOLD MRI in non-obese type 2 diabetic rats. Physiol Rep. 2016;4:e12890. doi: 10.14814/phy2.12890. [DOI] [PMC free article] [PubMed] [Google Scholar]

[R32] 32.Morikawa S, Kido C, Inubushi T. Observation of rat hind limb skeletal muscle during arterial occlusion and reperfusion by 31P MRS and 1H MRI. Magn Reason Imaging. 1991;9:269–274. doi: 10.1016/0730-725x(91)90411-e. [DOI] [PubMed] [Google Scholar]

[R33] 33.Morikawa S, Inubushi T, Kito C. Heterogeneous metabolic changes in the calf muscle of the rat during ischemia-reperfusion: in vivo analysis by 31P nuclear magnetic resonance chemical shift imaging and 1H magnetic resonance imaging. Cardiovasc Surg. 1993;1:337–342. [PubMed] [Google Scholar]

[R34] 34.Guruswami V, Sinop AK. Proceedings of the ACM-SIAM Symposium on Discrete Algorithms. Kyoto, Japan: 2012. Optimal column-based low-rank matrix reconstruction; pp. 1207–1214. [Google Scholar]

[R35] 35.Ma C, Lam F, Liu Q, Liang ZP. Proceedings of the ISMRM. Singapore: 2016. Accelerated high-resolution multidimensional 1H-MRSI using low-rank tensors; p. 379. [Google Scholar]

[R36] 36.Christodoulou AG, Shaw JL, Sharif B, Li D. Proceedings of the ISMRM. Singapore: 2016. A general low-rank tensor framework for high-dimensional cardiac imaging: Application to time-resolved T1 mapping; p. 867. [Google Scholar]

[R37] 37.Christodoulou AG, Redler G, Clifford B, Liang ZP, Halpern HJ, Epel B. Fast dynamic electron paramagnetic resonance (EPR) oxygen imaging using low-rank tensors. J Magn Reson. 2016;270:176–182. doi: 10.1016/j.jmr.2016.07.006. [DOI] [PMC free article] [PubMed] [Google Scholar]

PERMALINK

High-Resolution Dynamic 31P-MRSI Using a Low-Rank Tensor Model

Chao Ma

Bryan Clifford

Yuchi Liu

Yuning Gu

Fan Lam

Xin Yu

Zhi-Pei Liang

Abstract

Purpose

Methods

Results

Conclusions

INTRODUCTION

THEORY

Signal Model

Data Acquisition

FIG. 1.

Image Reconstruction

Subspace Estimation

Reconstruction

METHODS

Static and Dynamic 31P-MRSI at 3T

Dynamic 31P-MRSI at 9.4T

RESULTS

FIG. 2.

FIG. 3.

FIG. 4.

FIG. 5.

FIG. 6.

FIG. 7.