Accelerated 4D Quantitative Single Point EPR Imaging Using Model-based Reconstruction

k-Space Sampling Strategy

An incoherent sampling trajectory is a necessary component for sparsity-promoting reconstruction techniques where noise-like aliasing artifacts are desirable. The trajectory used herein utilizes conjugate symmetry and randomized sampling. Owing to pure phase encoding and low B0 field (10 mT), image phase in SP-EPRI is more tractable than in MRI (12). The reconstructed images in Figure 1 show a high degree of conjugate symmetry is possible in SP-EPRI. Therefore, we exploited conjugate symmetry when prescribing randomized sampling patterns by avoiding all symmetric points. Further, we designed a hierarchical random sampling scheme to effectively maximize the number of k-space samples in an undersampled acquisition.

Figure 2-a shows an example of the proposed hierarchical random sampling scheme, where k-space is segmented into 3 zones. Zone 1 is a fully sampled central region. Zone 2 is a region where k-space is undersampled by a factor of 2, which is converted to a fully sampled region after the application of conjugate symmetry. This large fully sampled region is required to maximize the performance of model-based reconstruction explained in later section. Zone 3 is more sparsely sampled than Zone 2. In Zone 2, k-space is uniform randomly sampled, whereas in Zone 3, Gaussian random sampling is applied. After sampling, the sampled points in Zone 2 and Zone 3 are used to estimate the samples at the symmetric position by the application of conjugate symmetry. In this study, 10% of matrix size was used for the width of Zone 1, and the size of Zone 2 was set to 50%, 40%, and 33% respectively for 61×61×61, 95×95×95, and 127×127×127 dataset. 40% of the matrix size was used for the standard deviation of the Gaussian random sampling in Zone 3.

Bilateral k-space extrapolation

Recently, we proposed a gridding technique for SP-EPRI, called k-space extrapolation (KSE), to achieve equal FOV reconstruction with time-invariant Gibbs ringing and thereby enable more reliable pixelwise T₂* estimation (6). Unfortunately, this method is not suitable for highly undersampled k-spaces since high frequency regions will be too sparsely sampled and there exist too few samples to be extrapolated from the later phase encoding time delays to earlier delays. Therefore, this method has been improved by performing k-space extrapolation bilaterally as shown in Figure 2-b.

When performing bilateral k-space extrapolation (bi-KSE) with a hierarchical sampling pattern, Zone 1 and Zone 2 do not need bilateral extrapolation since those regions are fully sampled (Zone 2 becomes equivalent to fully sampled after the application of conjugate symmetry). Therefore, Zone 1 and Zone 2 are extrapolated using unilateral/backwards KSE, while Zone 3 is extrapolated bilaterally as described in Figure 2-b. The extrapolated k-space samples need to be scaled to harmonize with the target k-space, compensating for T₂* decay. This scaling factor can be approximated by simply referring to the center of k-space. However, when performing k-space extrapolation within the central region (low spatial frequencies), small errors in the estimated scaling factors can lead to severe distortion in image reconstruction and parameter estimation. To address this issue, the scaling factor is refined iteratively so that the Gibbs ringing pattern of the extrapolated image becomes similar with the reference image (the one with the highest bandwidth). To evaluate the similarity between two images, the numerical gradient of the images were used where gradient values near strong edges were discounted to purely consider changes resulting from Gibbs ringing. The numerical gradient images were obtained using the central differencing scheme (23). The scaling factor was adjusted by evaluating the error function that is the l₁-norm of error between the target and reference gradient image, using the Nelder-Mead simplex search algorithm as implemented by MATLAB (The Mathworks, Inc., Natick, MA).

After bilateral k-space extrapolation, equal FOV images are reconstructed by applying convolution gridding (24,25). Iterative density correction was used before gridding (26), which is an indispensable process for reliable quantification, because spatial sampling density changes at each reconstructed time delay. As the reconstructed images show in Figure 3-a, bi-KSE dramatically improves the quality of images reconstructed with undersampled k-spaces, especially for images at later time delays (1300 ns) where reconstruction error is greatly reduced (Figure 3-b). The proposed bi-KSE method dramatically increases available k-space samples (74,379 vs 1,191 at 1300 ns for bilateral and unilateral k-space extrapolation, respectively) and provides improved image quality compared to previous techniques. However, not all reconstruction error is removed. Therefore, we employed PCA-based reconstruction to exploit the rich spectral domain of the SP-EPRI dataset to further improve image quality and resulting T₂* estimation, as explained below.

Model-based reconstruction

To evaluate model based reconstruction in highly undersampled SP-EPRI, we have utilized PCA-constrained reconstruction (19). In PCA-based reconstruction, a training matrix whose columns consist of training data, which consist of time series for all possible FID signals, are used to obtain principal components. The eigenvectors or singular vectors obtained from the training matrix (by using eigenvalue decomposition or singular value decomposition) will span the signal subspace where true FID signals exist.

The training data is generated within a predefined T₂* range and is a factor that affects the performance of PCA-based reconstruction. Utilizing a T₂* range that does not encompass the T₂*’s of targeted objects may result in a false signal space. In hypoxia imaging EPRI, the T₂*’s of Oxo-63 are distributed within a small range, approximately 400-650 ns corresponding to a dissolved oxygen level between 5% and 0%, respectively. To encompass all possible oxygen levels (0-100%), a T₂* range of [1ns, 700ns] was used in this study.

With the estimated T₂* range [p,q] (p<q), a training matrix can be composed using the independent monoexponential curves as Equation 1.

Then, D is decomposed by SVD to yield singular vectors, and L significant singular vectors (L=3 or 4 was used in this study) are selected to compose matrix B̂. PC coefficient matrix M can be obtained by Equation 2, whose n-th column vector represents PC coefficient map corresponding to n-th PC.

, where I is an image matrix containing vectorized initial images, T is the number of images, and N is the length of image vector. Now, the model-based reconstruction problem can be formulated by the following equation:

, where FT_j denotes the operator for discrete Fourier transform with image at time delay j, $\overrightarrow{K_J}$ is measured k-space vector at time delay j, P_i(M) is penalty based on discrete wavelet transform (Daubechies 4, Wavelet Toolbox in MATLAB 2011b) and total variation (27) of PC coefficient maps, and λ_i is the Lagrange multiplier that is selected differently with each PC coefficient map. The penalty term is calculated as in Equation 4:

, where ψ is a DWT transform, TV is total variation, M_i is i-th PC coefficient map, and α is tuning constant between two objectives.

In Equation 3, the first and the second summation term in the brace represent data fidelity and penalties respectively. In practice, calculation of the data consistency term requires significant computational power due to repeated gridding and inverse gridding steps (to deal with non-Cartesian k-space samples accumulated by bi-KSE), especially when dealing with large numbers of samples in 3D imaging. To improve the speed of reconstruction, we enforced data fidelity in the image domain rather than in the k-space domain as the following equation shows:

, where (MB̂^T)_j and $\overrightarrow{I_J}$ denote j-th column vector in respectively updated and initial image matrix. This approximation is possible owing to the hierarchical random sampling pattern which provides good quality image reconstruction after the application of bilateral k-space extrapolation (as seen in Figure 3 and Figure 6-g). In experiments not presented herein, image-based consistency performed as well as conventional k-space methods, likely due to the strong performance of bi-KSE alone for undersampled acquisitions.

This optimization problem was solved using the nonlinear Polak-Ribiere Conjugate Gradient algorithm initiated with steepest gradient and backtracking line search with contraction factor 0.4. Once the PC coefficient maps are optimized, the image sequence can be reconstructed by performing a linear combination of PC coefficient maps and PC vectors as Equation 6.

In this method, the image sequences were normalized by the maximum value of the initial image in the sequence. By doing this, we were able to limit the effective range for parameter setting (α, λ₁, λ₂, λ₃, and λ₄), which can be generally applied in different experiments without dependence on any scaling factors between datasets. The effective parameter ranges 0.1 ≤ α ≤ 1.0, and 0.1S_k ≤ λ_k ≤ 0.5S_k were empirically found and used, where ${\left(S_k\right)}_{k=1}^4=(1,2,5,10)$ for the k-th PC coefficient map. The selected λ_k is then scaled by ∥M_k∥/∥M₁∥ to compensate for the scale difference between PC coefficient maps. Note that larger λ_k is used for the less significant PCs that commonly contain more noisy data. By using larger λ_k we can promote sparsity and smoothness in the corresponding PC coefficient map and thereby suppress noise more effectively.

Simulation

To evaluate the capability of the proposed method for T₂* estimation with undersampled k-space, a computer simulation and a phantom experiment using SP-EPRI were performed. For the computer simulation using MATLAB (The Mathworks, Natick, MA), synthesized SP-EPRI images with much higher resolution (187×187×187) than conventional acquisitions (e.g., 19×19×19) were generated based on the 3D T₂* map shown in Figure 7-a. From the 3D T₂* map, FID curves were simulated from 1ns to 1800ns using a sampling rate of 5ns and a dead time of 300ns. To simulate SPI encoding, inverse gridding was used to sample k-space with a 61×61×61, 95×95×95 or 127×127×127 matrix with a spreading Dirac comb function to simulate the zoom-in effect due to constant gradients. For 2D phantom experiments, only the central slice of the phantom was used.

Simulation - Undersampling vs. Reduced Averaging

The voxel size for the proposed acquisitions is 25~220× smaller than previous SP-EPRI acquisitions, which results in a large reduction in SNR. Even for conventional resolutions, a large number of averages (1000~8000) is typically applied in SP-EPRI to improve SNR. Therefore, it may be also possible to accelerate imaging by simply reducing the number of averages rather than undersampling k-space; however, this is not expected to perform well due to the already low SNR of the data. To verify this, a simulation was performed to compare the proposed undersampled k-space acquisition to reduced average data, with and without the proposed PCA-constrained reconstruction. The reduced average method was simulated by adjusting noise levels to decrease SNR by a factor of $\sqrt{R}$ (R: acceleration factor attained by reducing average). A 127×127 2D digital phantom was generated as explained above to simulate comparable undersampled and reduced averaging datasets with R=8.

The 2D simulation was performed at the critical SNR limit for the proposed method, where critical SNR is the lower SNR limit in the region of shortest T₂* at time delay 1300 ns that still enables reasonable parameter estimation. For these experimental conditions, this was determined to be approximately 3, as shown in Figure 4. Simulated data with an SNR of 3 was undersampled with R=8 and processed by the proposed method explained in the above sections. For the reduced average method with R=8, the SNR was reduced by a factor of $\sqrt{8}(\mathrm{SNR}=1.06)$, and then the data were fully-sampled. The sampled k-spaces were processed with the conventional unilateral KSE and PCA-based reconstruction. 101 consecutive k-spaces from 800ns to 1300ns were used for reconstruction.

Simulation - Resolution and Acceleration

The fully sampled k-spaces were retrospectively undersampled with different acceleration factors (R) by using the hierarchical random sampling explained in the above section. The data were generated with low SNRs to highlight the performance of the proposed method (5 for 127×127×127, 7 for 95×95×95, and 9 for 61×61×61 k-spaces, when measured at 1100 ns). Data were acquired at R=1, 4, 6, 8 for 61×61×61 k-space; R=8, 12, 15 for 95×95×95; R=15, 30, 60 for 127×127×127 k-space.

Phantom experiment

Details regarding the specifications of the EPRI spectrometer have been previously published (5,28). For the tube phantom experiment, 3D EPR data was obtained using three-tube phantom comprised of three tubes containing 2 mM Oxo-63 (GE Healthcare, Waukesha, WI) bubbled with 0%, 2% and 5% oxygen, respectively. Data was encoded using three orthogonal phase-encoding gradients incrementally ramping in 61 equal steps from −40 to 40 mT m^-1, resulting in 61×61×61 phase-encoding steps. 581 data points were encoded with a sampling period of 5ns after the minimum RF recovery dead time (530 ns). 8000 averages per phase encoding point and an interpulse delay (TR) of 10 μs was used. A total of 28,373 points (approximately 8-fold undersampling) with a 61×61×61 matrix size were phase-encoded using the hierarchical random sampling strategy.

Data Processing

Bi-KSE was performed in the peripheral region (Zone 3), and unilateral/backwards KSE was independently performed in the central region (Zone 1 and Zone 2) in a reversely cascading manner, as depicted in Figure 5-a. The extrapolated k-spaces from each phase encoding time delay were individually gridded into Cartesian k-space to generate equal FOV images (Figure 5-b). Then, these images were used as input to the PCA-based reconstruction (Figure 5-c). With a sequence of final images reconstructed by using above-explained methods, T₂* was fit using a traditional T₂* relaxation model for the magnitude of pixelwise transverse magnetization (Figure 5-d). Linear least-square curve fitting was applied to the log-linearized FID curve in which the slope and the y-intercept represent respectively -1/T₂* and log(M₀).

Simulation Results – Undersampling vs. Reduced Averaging

The reconstructed images and the resultant T₂* maps obtained using reduced average or undersampling are shown in Figure 6. Reconstructions with gridding are shown in 6a, 6b, 6c for reduced average data, and 6f for undersampled data. KSE is applied to the reduced average data in 6d and 6e, and bi-KSE is applied to the undersampled data in 6g and 6h. PCA-based reconstruction is applied to the reduced average data in 6c and 6e, and to the undersampled data in 6h. As seen in the reconstructed images and the resulting estimated T₂* maps, undersampling (6h) outperforms reduced averaging (6c and 6e).

Simulation Results – Resolution and Acceleration

To evaluate how accurately the proposed method estimates T₂* with undersampled k-space data, we performed computer simulations using synthesized 3D data consisting of various T₂*s (Figure 7-a). In addition to the proposed method, a conventional k-space extrapolation method with full sampling (R=1) was also implemented for comparison. 81 consecutive k-spaces from 700 ns to 1100 ns with a 5 ns time interval were used for reconstruction. Figure 7-b shows the estimated T₂* maps obtained using undersampled k-spaces with 61×61×61 gradient steps with acceleration factors of R=1, 4, 6, and 8. Table 1 shows the corresponding quantitative result. Our method enabled accurate T₂* estimation up to 8-fold acceleration at this matrix size. However, the small segment E tends to be blurred due to the large voxel size. Figure 8 depicts how the proposed method works with larger matrix sizes. As seen, when higher matrix size is used, the overall performance of our method tends to be improved at the same acceleration factor (R=8) due to the larger central region and higher expected transform sparsity. Moreover, good quality T₂* estimation with even higher acceleration factors (R=15 for a 95×95×95 k-space, and R=30 for a 127×127×127 k-space) is attainable.

Phantom Experiment Results

Figure 9 shows the T₂* map estimated from prospectively undersampled k-spaces (R=8). 91 consecutive k-spaces from 1530 ns to 1980 ns with a 5 ns time interval were used for reconstruction. As seen in Figure 9, the proposed method enables reliable parameter estimation with reduced sampling. The estimate was 684.45 ± 39.10 ns, 659.79 ± 31.11 ns, and 591.10 ± 25.52 ns for tube bubbled with 0%, 2%, and 5% oxygen, respectively, which lie within the our expected range. The estimated numbers show 4%-6% standard deviation, which is presumably due to the system noise in the current EPR scanner. Fitted %oxygen-R₂* curve showed slope of 4.68×10^-5 (ns^-1 %oxygen^-1) and y-intercept of 1.45×10^-3 (ns^-1) with R²=0.9662.