Fast Padé Transform Accelerated CSI for Hyperpolarized MRS.

The fast Padé transform (FPT) is a method of spectral analysis that can be used to reconstruct nuclear magnetic resonance spectra from truncated free induction decay signals with superior robustness and spectral resolution compared with conventional Fourier analysis. The aim of this study is to show the utility of FPT in reducing of the scan time required for hyperpolarized 13C chemical shift imaging (CSI) without sacrificing the ability to resolve a full spectrum. Simulations, phantom, and in vivo hyperpolarized [1-13C] pyruvate CSI data were processed with FPT and compared with conventional analysis methods. FPT shows improved stability and spectral resolution on truncated data compared with the fast Fourier transform and shows results that are comparable to those of the model-based fitting methods, enabling a reduction in the needed acquisition time in 13C CSI experiments. Using FPT can reduce the readout length in the spectral dimension by 2-6 times in 13C CSI compared with conventional Fourier analysis without sacrificing the spectral resolution. This increased speed is crucial for 13C CSI because T1 relaxation considerably limits the available scan time. In addition, FPT can also yield direct quantification of metabolite concentration without the additional peak analysis required in conventional Fourier analysis.

Introduction

Hyperpolarized Magnetic Resonance (MR) with Dissolution-Dynamic Nuclear Polarization is a clinically emerging technique, and it shows great promise in investigating diseases of metabolic dysregulation such as cancer and heart diseases (1). Chemical shift imaging (CSI) is a well-known and powerful tool for the investigation of metabolic processes following the injection of hyperpolarized MR media. Because the T1 relaxation rate of hyperpolarized 13C places an ultimate limit on the time of acquisition of CSI data (forming a window of 50-80 seconds in vivo) (2, 3), spectroscopic 13C imaging is typically performed with a multipoint Dixon approach (4) to limit the number of time points required for spectral analysis. However, this approach requires a priori knowledge of the number and spectral location of spectral peaks expected to be present in the spectrum. If other signals are present at unexpected resonant frequencies (ie, unexpected reactions occur), these multipoint methods are prone to producing incorrect results. In contrast, full-spectrum CSI yields a complete MR spectrum for each image voxel, with zero to few prior assumptions. However, the acquisition time necessary to obtain full spectra with adequate resolution from spatially resolved voxels represents a significant limitation in hyperpolarized 13C imaging.

Numerous novel methods have been proposed for gaining higher spatial and temporal resolution at the expense of spectral resolution (4–10). Undersampling the spectral dimension can considerably reduce the amount of data that is required to be acquired and hence the acquisition time; this is typically offset by the introduction of uncertainties in the spectral domain of the acquired data, which often require prior spectral knowledge to reconstruct, causing a potential for misinterpretation. Consequently, despite an increasing number of alternatives, the CSI sequence is still considered a gold standard sequence for hyperpolarized imaging experiments and a preferred starting point for testing experimental hyperpolarized setups (7, 11–18).

The CSI sequence makes few prior assumptions about obtained spectral frequencies and is therefore both more robust experimentally and less prone to misinterpretation than faster alternatives. This flexibility in the acquisition and processing makes the CSI sequence a valuable tool for clinical and preclinical investigations. The per-shot acquisition time for optimal signal-to-noise ratio (SNR) is typically taken to be in the range of 1 to 2 times the T2* value, maximizing acquisition of the free induction decay (FID) and minimizing the sampling of background noise. Owing to the need of fully sampling an individual FID for each resolved voxel, CSI is traditionally limited in temporal resolution. Reducing the imaging time necessary to acquire full spectra with adequate spectral resolution will increase the temporal resolution of the CSI experiment.

After acquisition of a CSI data set, the identity and relative concentrations of hyperpolarized metabolites are often determined by measuring the individual peak heights, integrating peaks, or fitting models to the resonances in the spectra (7, 13, 14, 19) obtained after Fourier transforming CSI data. These methods are robust when peaks are well separated and nonoverlapping. Hankel singular value decomposition methods have also been proposed in various papers (20). Here, we present an alternative approach, using the fast Padé Transform (FPT) as opposed to the fast Fourier transform (FFT), to directly quantify both simulated CSI data and data obtained from the spectral imaging of hyperpolarized [1-13C]pyruvate in the rat kidney.

Theory

FFT is the traditional method for transforming data in the time domain to the frequency domain. If the Nyquist criterion is satisfied, it is a robust and efficient tool to obtain the spectral envelopes of metabolic data sets when no or very little prior knowledge is present. However, the interpretation and quantification of the results are limited when the FID is truncated because of noise or rapid acquisition, as the artificial zero-filling of the FID results in a sinc-shaped interpolation with oscillating side-lobes around each resonance present. For any pair of closely spaced resonances, truncating the acquisition of the FID can produce spurious peaks or dips in the resulting spectrum. Correspondingly, parametric approaches are often taken to minimize the effect of this interpolation and quantify the spectra by fitting Lorentzian or Gaussian lineshapes to their envelopes in a least-squares fashion. Popular spectral and time-domain fitting algorithms such as AMARES (21), LCModel (22), and VAPRO (19) use user-supplied prior information to find the spectral parameters of interest to the Fourier transformed data fitting, either a predicted FID in the time domain or the envelope of the spectrum in the frequency domain. Such fits are not unique, but are robust and straightforward to interpret (23).

An alternative approach to fitting spectral data is the direct extraction of frequency and amplitude information from the FID by methods other than FFT by the construction of objects called parametric estimators, and their utility in MR is well known. In particular, FPT has been shown to enhance the resolution of nuclear magnetic resonance (NMR) spectra and to be less affected by truncation artifacts and baseline modulations than FFT alone (23–26). The derivation of FPT is further described.

Let be a discretely sampled signal, where c is the (n + 1)th value of the FID composed of N data points. Then, in general, the discrete Fourier components of c is given as follows: where τ is the sampling delay between acquired points. Defining z to z := exp(iωτ), we arrive at the general frequency spectrum defined as follows:

The Padé transform constructs a Padé approximant to S (z) before exploiting analytical properties of rational polynomials to extract quantities of interest such as amplitudes and frequencies from this polynomial. The Padé approximant of S is defined as a quotient of polynomials, A (z)/B (z), such that for as many k = 0, 1, … as possible, there holds the following: where are complex polynomials of degree P and Q, respectively, and where we set P = Q = N/2. The set of unique coefficients {a} can be found by noting that A (z) = B (z) S (z) and equating equal powers in z, yielding 2 systems of linear equations. Singular value decomposition yields {bj} efficiently, and {aj} can then be explicitly computed given the data. For the systems considered in hyperpolarized MR experiments, where N typically is 1024-4096, FPT reconstruction does not present a significant computational challenge.

Once A (z) and B (z) are known, it can be shown that the k roots of B (z) are related to the resonances in the spectrum as follows: for which several efficient numerical algorithms exist. The corresponding complex amplitudes of each ωk are derived from the Cauchy residue of the rational approximant A (z)/B (z) evaluated at the given root k, which is known analytically as follows:

Accordingly, the amplitude and spectral position of resonances can be recovered, and the spectrum is generated in any desired mode (absorption, dispersion, and magnitude) by using Heaviside partial fractions (27), that is, by the expansion.

It can be shown that Re (ωk) and Im (ωk) are the position and width of the resonance k, and dk and arg (dk) are its amplitude and phase, respectively. Thus, FPT can find every peak parameter of interest for a resonance, without ever using the Fourier spectrum or a fitting procedure. In the case that B (z) contains unique roots, the returned spectrum consists of several Lorentzians; if B (z) has repeating roots, then the resulting spectrum contains spectrally overlapping resonances with a non-Lorentzian lineshape. Depending on the shim of the region of interest (ROI), these multiple resonances may be spurious. In addition, the Padé approximation has the advantageous property that the analytical approximant constructed can be used for spectral extrapolation inherently, analogously to zero-filling the FID before calculating FFT (28).

Experimentally acquired spectra contain noise, whose presence can manifest in FPT formalism as common zeros of A (z) and B (z), is referred to as Froissart doublets for historical reasons, which can yield spurious signals (ie, artifacts) in the approximation. Froissart doublets are poles and zeroes with a numerical distance of zero in the noiseless setting, and in real experimental data, close to zero (23, 24, 26, 29) with a large residue (30). The doublets can also be observed through the iteration of FPT as nonconverged frequencies and can be removed by a threshold in relation to the background noise levels. Here, the noise threshold contribution was set as amplitudes occurring within 5 standard deviations of the noise as manually identified in the edges of the spectrum (ie, the outer 10% of the bandwidth). An iterative approach, where the number of sampled points K is varied with K = Q, can identify spurious signals, as noisier signals have a higher likelihood of changing parameters in-between iterations (23, 24, 26).

Methods

Simulated Data

Simulated time-domain data from an NMR spectrum with 5 peaks typical of a [1-13C] pyruvate experiment was processed with FFT with 256 points over a 4000 Hz bandwidth and a truncated FPT with 16 points (K = Q = 8) and the same bandwidth, interpolated to 256 points. Frequency components found by FPT was reconstructed both with a Gaussian and a Lorentzian model for comparison.

Phantom Data

[1-13C]urea and [1-13C]acetate phantoms (with a spectral separation of 1890 Hz at 9.4 T) were scanned to explore image properties of the accelerated acquisition using a 9.4 T horizontal bore magnet (Agilent, Yarnton, UK) with VnmrJ 4.0A (Agilent, Santa Clara, California). A dual-tuned 13C/1H volume rat coil (Doty Scientific, Columbia, South Carolina) was used for 1H and 13C magnetic resonance imaging. A standard slice-selective 2-dimensional 13C CSI sequence was used for imaging the phantoms. Parameters were as follows: flip angle = 10°, a Cartesian k-space trajectory, matrix = 32 × 32, repetition time/echo time = 200 ms/0.67 ms, field of view = 60 × 60 mm2, average = 4, spectral width = 6000 Hz, number of points = [256 128 64 42 32 16] complex points, and an axial slice thickness of 20 mm. FPT reconstruction of the accelerated data was compared with the iterative decomposition of water and fat with echo asymmetry and least-squares estimation (IDEAL) method to assess the image quality (31).

Magnetic Resonance Spectroscopy Data

Experimental magnetic resonance spectroscopy (MRS) (n = 6) data [obtained from an earlier study (14)] were collected from rat kidneys using a 4.7 T horizontal bore magnet (Oxford Instruments, Oxford, UK) equipped with a Varian Direct Drive console and VnmrJ 2.3A (Agilent, Santa Clara). A bicarbonate phantom was placed next to the rat. A dual-tuned 1H/13C volume coil with a 13C 4-channel phased array receive coil (Rapid Biomedical, Würzburg, Germany) was used for 1H and 13C magnetic resonance imaging and MRS, respectively. The kidneys were localized by a standard gradient-echo sequence, and a slice covering both kidneys was manually shimmed. A standard slice-selective 2-dimensional 13C CSI sequence was used for hyperpolarized [1-13C] pyruvate imaging. The parameters were as follows: flip angle = 10°, a full centric circular k-space trajectory, matrix = 16 × 16, repetition time/echo time = 75 ms/0.65 ms, field of view = 60 × 60 mm2, spectral width = 4000 Hz, number of points = 256 complex points, and an axial slice thickness of 20 mm, covering both kidneys. The script for FPT CSI processing was implemented in MATLAB (MathWorks, Natick, Massachusetts) available at http://www.mr.au.dk/software. Both truncated FPT and standard FID-CSI data were processed in MATLAB, and spatial dimensions were apodized with a Hamming function and zero-filled to a 32 × 32 matrix. Different levels of truncation of the FID data were investigated using the initial 1/2, 1/4, 1/6, 1/8, and 1/12 of the total data. The spectral domain was processed as sum-of-squares of the 4 coil channels. The metabolite maps were imported to OsiriX, and an ROI analysis was performed. The metabolite maps were normalized relative to the sum of all metabolites. Statistical plots were produced in GraphPad Prism (GraphPad Software, Inc. La Jolla, San Diego, California).

Results

FPT of the noiseless simulated data identifies genuine signals (Figure 1) even from extremely truncated signals. No prior knowledge is required for identification of the genuine frequencies, as the Froissart doublet assumption is valid in the noiseless situation. The Lorentzian model produces better results as a model for the extracted metabolites when compared with the Gaussian model (Figure 1B). One Froissart doublet is seen as a Q and P root at the same position (Figure 1C).

Figure 1.

Simulated [1-13C] pyruvate signal, processed with fast Fourier transform (FFT) with 256 points over 4000 Hz and a truncated fast Padé transform (FPT) by using the first 16 points (K = Q = 8), zero-filled to 256 points over 4000 Hz, showing similar overall resolution (A). The individual components of the FPT show a good reasonable fit with K = Q − 1 = 7 frequencies with the Lorentzian model. Using a Gaussian model, the spectrum show less precise results (B). The Froissart doublets show 2 spurious signals inside the unit circle, where the zk,Q = zk,P, ie, dk = 0 (C).

In Figure 2, the comparison between accelerated data reconstructed by FPT and IDEAL can be seen. FPT reconstructs similar images compared with the IDEAL method, and it is able to separate the 2 phantoms even with very high acceleration. Phantom signal and noise levels are consistent in FPT reconstruction, and shapes features are preserved.

Figure 2.

[1-13C]urea and [1-13C]acetate phantoms (with a spectral separation of 1890 Hz at 9.4 T) reconstructed with FPT and iterative decomposition of water and fat with echo asymmetry and least-squares estimation (IDEAL) are shown in the images on a logarithmic scale. The images with progressively reduced numbers of points acquired show a great similarity between the 2 methods.

Truncating the in vivo FID caused some inaccuracy in the identification of frequencies' peaks by FPT (Figure 3). As shown in Figure 3, FPT was more resilient to FID truncation. The apparent degradation of the NMR spectrum caused by FID truncation was much less with FPT analysis, compared with FFT analysis, with spectral lineshapes notably preserved by FPT and not FFT. In addition, quantification of metabolite concentrations in the frequency domain with conventional FFT analysis requires an additional computation to be performed on the NMR spectrum (eg, peak fitting or peak integration).

Figure 3.

Representative data set comparing FFT and FPT performance with varying amounts of truncated data. Spectrum from region of interest (ROI) position is shown with truncation 1/2, 1/4, 1/6, 1/8, and 1/16 of the original data.

The ROI analysis of the metabolite maps formed by FPT was consistent with maps created by full-dataset FFT as the degree of signal truncation increases (Figure 4). However, the ROI analysis of the MRS data showed that increasing truncation of the FID increases the contribution of signal from pyruvate in the kidneys by both FPT and FFT. On FFT maps, the increased FID truncation resulted in increasing errors in concentration estimation of the different metabolites. In comparison, FPT maps were relatively unaffected by FID truncation. To assess stability against truncation, an ROI was placed in the kidney and in the urea phantom. The difference in bicarbonate and urea signal between the 2 sites was then measured. As seen in the urea phantom signal (Figure 4), FFT urea signal decreases with truncation, whereas FPT shows signal stability even with high truncation factors. Figure 5 shows a significant (P < .05) difference between FFT and FPT in the bicarbonate signal from this truncation at (1/4, 1/6, 1/8, and 1/12) parts of the original signal. In the case of urea signal, FFT and FPT showed a significant difference from truncation at (1/6, 1/8, and 1/12) parts of the original signal (Table 1). The filtering features of FPT are initiated in the fully sampled data set, which can be seen in Table 1, and show a slightly higher SNR compared with FFT.

Figure 4.

ROI analysis was obtained from 3 locations (kidney, background, and urea phantom) in both FFT and FPT images with increasing truncation, 1/2, 1/4, 1/6, 1/8, and 1/16 parts of the full dataset (x-axis, 1-6). The FPT shows an increasing signal separation for the 3 selected metabolites, and FFT shows a decrease and has severe truncation error on the bicarbonate.

Figure 5.

Difference between kidney and phantom mean bicarbonate ROI signal in the bicarbonate image (visual representation of Table 1) (A). Difference between kidney and phantom mean ROI signal in the urea image (B). As the degree of data truncation increases, signals from the different metabolites are blurred by FFT analysis. This effect is significantly diminished with FPT analysis (P < .05).

Table 1.

The Difference Between FFT and FPT in Urea and Bicarbonate

	Number of Samples	256	128	64	42	32	21
Bicarbonate	P value	0.556	0.566	0.046	0.027	0.044	0.000
	Difference	−0.008	−0.008	−0.033	−0.040	−0.037	−0.106
	SE of difference	0.012	0.013	0.015	0.015	0.016	0.017
	FFT: Mean ± SD	0.020 ± 0.017	0.008 ± 0.025	−0.011 ± 0.033	−0.019 ± 0.032	−0.015 ± 0.035	−0.090 ± 0.037
	FPT: Mean ± SD	0.028 ± 0.025	0.015 ± 0.021	0.022 ± 0.014	0.021 ± 0.020	0.022 ± 0.018	0.016 ± 0.021
Urea	P value	0.130	0.312	0.041	0.005	0.002	0.000
	Difference	−0.076	−0.057	−0.123	−0.203	−0.219	−0.335
	SE of difference	0.046	0.053	0.053	0.056	0.053	0.045
	FFT: Mean ± SD	0.272 ± 0.087	0.286 ± 0.100	0.263 ± 0.095	0.207 ± 0.085	0.183 ± 0.076	0.117 ± 0.066
	FPT: Mean ± SD	0.348 ± 0.072	0.342 ± 0.084	0.387 ± 0.088	0.411 ± 0.107	0.402 ± 0.105	0.452 ± 0.087

Notes: Mean ROI signal in urea and bicarbonate detected using FFT and FPT methods.

Discussion

FPT offers a novel method for reducing the acquisition time for spectroscopy by a factor of 2-6, without decreasing the apparent resolution of the NMR spectrum, thereby significantly accelerating 13C CSI. A 16-second total scan time can be reduced to 3 seconds by reducing FID by 6, as the sampling time is a major determinant of the minimum repetition time. A reduction to 50% sampling reduces the scan time to 8 seconds. FPT can be an alternative to subsampled FFT methods, which require prior knowledge of metabolites and produce inaccurate results when unanticipated signals are present. There are benefits to obtaining full-spectral data compared with specific metabolic imaging with, for example, a spectral–spatial pulse as follows: spectroscopic approaches are relatively immune to bulk frequency offsets or B0 inhomogeneity, whereas either of these effects can prevent the successful acquisition of a spectral–spatial dataset. Using FPT allows “full spectrum truncated” 13C CSI, which is not possible with conventional FFT analysis. FPT assumes a Lorentzian lineshape for each signal in the frequency domain, resulting in nonLorentzian signals being represented by several components, a biophysically meaningless fit; however, a simple integral over a specified region incorporated these components. The results from the experimental MRS data give a good indication of the stability provided by FPT. Unwanted bleeding of frequency signal to nearby metabolites occurs in FFT analysis of truncated data but not in FPT, which, in the future, may be a key feature in choosing FPT. The model-based IDEAL reconstruction approach shows superior acceleration capabilities compared with FPT; however, this is at the expense of little flexibility post acquisition and the necessity for a prior knowledge. The difference in SNR between IDEAL and FPT originates from the difference in integration of the peak (FPT) and fitting the exact peak (IDEAL).

A priori and a posteriori knowledge is easily incorporated into FPT during the sorting of the Froissart doublets, allowing simple and intuitive filtering of the spurious spectral components. An iterative filtering process can improve the SNR and the genuine frequencies' recognition by varying K, using the Froissart doublets and the robustness in the genuine peaks, with regard to varying the frequency, damping, and area, thereby filtering out frequencies that are not present in all iterations (29, 32). The method is similar to the well-known Hankel singular value decomposition methods. However, the unique utilization of the Froissart criteria makes FPT approach an interesting alternative to previous methods, with the potential of reducing the overall experimental time, without losing the flexibility in the acquisition and processing.

The main advantages of FPT are the direct quantification of spectral data and the possibility of considerably decreasing acquisition time. Future directions for research based on this work include implementation of a no prior knowledge of metabolic selection filter for 13C spectroscopy and incorporation of compressed sampling methods to further decrease the required acquisition time and combine FPT with acceleration methods such as parallel imaging (12) and compressed sensing (33), allowing further reduction in scan times.

Conclusion

We have implemented FPT as a method for estimating metabolite concentrations in hyperpolarized 13C chemical shift imaging. The short T1 of hyperpolarized 13C severely limits the time available for signal acquisition, making conventional CSI difficult or impossible. However, with FPT, spectra can be reconstructed from truncated FID signals without the loss of spectral resolution. This allows shortening of the data acquisition time, making 13C CSI more feasible. We have shown, with numerical simulations and in vivo 13C CSI data, that FPT gives more accurate quantification of metabolites with truncated data acquisition than conventional Fourier analysis.