Compressed Sensing for fMRI: Feasibility Study on the Acceleration of Non-EPI fMRI at 9.4T.

Conventional functional magnetic resonance imaging (fMRI) technique known as gradient-recalled echo (GRE) echo-planar imaging (EPI) is sensitive to image distortion and degradation caused by local magnetic field inhomogeneity at high magnetic fields. Non-EPI sequences such as spoiled gradient echo and balanced steady-state free precession (bSSFP) have been proposed as an alternative high-resolution fMRI technique; however, the temporal resolution of these sequences is lower than the typically used GRE-EPI fMRI. One potential approach to improve the temporal resolution is to use compressed sensing (CS). In this study, we tested the feasibility of k-t FOCUSS--one of the high performance CS algorithms for dynamic MRI--for non-EPI fMRI at 9.4 T using the model of rat somatosensory stimulation. To optimize the performance of CS reconstruction, different sampling patterns and k-t FOCUSS variations were investigated. Experimental results show that an optimized k-t FOCUSS algorithm with acceleration by a factor of 4 works well for non-EPI fMRI at high field under various statistical criteria, which confirms that a combination of CS and a non-EPI sequence may be a good solution for high-resolution fMRI at high fields.

1. Introduction

Functional magnetic resonance imaging (fMRI) has had a wide impact in both the research and clinical community since its development. In conventional fMRI studies, positive blood oxygen level-dependent (BOLD) response signal is used as a measure to map neural activity in the brain [1], and the most common MR pulse sequence for acquiring BOLD fMRI images has been gradient-recalled echo (GRE) echo-planar imaging (EPI) due to its fast acquisition speed and high sensitivity to BOLD effect. However, this technique is susceptible to local magnetic field inhomogeneity and becomes sensitive to image distortion and degradation especially at high magnetic fields. Non-EPI sequences such as spoiled gradient echo or balanced steady-state free precession (bSSFP) can be used as an alternative tool for fMRI [2-9]; however, the major drawback of using these sequences for fMRI studies is the low temporal resolution compared to the typically used GRE-EPI.

One solution to overcome the low temporal resolution of non-EPI sequences is to adopt parallel imaging technique [10-12]. Although proven useful, the usage of parallel imaging results in reduced signal-to-noise ratio (SNR) due to the acceleration factor, the geometric factor of the different coil elements, and the k-space filling trajectory. The other solution to improve the temporal resolution of non-EPI sequences is to use compressed sensing (CS) [13-15]. CS theory states that it is possible to reconstruct an aliasing-free image even at sampling rates dramatically lower than the Nyquist sampling limit, as long as the nonzero signal is sparse and sampled incoherently. These requirements can be well satisfied in dynamic MRI, since arbitrary trajectories can be implemented to incoherently sample data and dynamic MR images can be sparsified due to high temporal redundancy [16, 17]. Recently, CS theory was successfully applied to dynamic MRI in a new algorithm called k-t FOCUSS by employing random sampling pattern in k-t space and by using various sparsifying temporal transforms such as Fourier transform (FT) and Karhunen-Loéve transform (KLT) to utilize the temporal redundancies [16, 18].

Though CS theory has gained attraction for its vast potential for MRI application, CS had been successfully applied to fMRI in only a few studies in the past. In most of the studies, CS was applied to GRE-EPI fMRI: ordinary GRE-EPI [19, 20] and spiral scan GRE-EPI [21]. Despite its application, GRE-EPI is generally known to suffer from the contribution of magnetic field inhomogeneity, which can degrade the performance of CS algorithms. Recently, it has been reported that application of CS to GRE-fMRI may increase statistical performance of activation detection [22]. Non-EPI sequences such as bSSFP or GRE may work better with CS algorithms than GRE-EPI, since the sequences utilize different RF excitations for each TR.

The signal dynamics of steady-state sequences such as bSSFP are known to be more complicated than other conventional sequences and may need careful examination prior to the application of CS to real data. Also, due to the large degree of freedom in CS application, it is important to understand the artifacts and effects related to CS reconstruction without any confounding factors from true physical artifacts. Thus, an extensive simulation study prior to actual CS application is required to reconstruct data appropriately, preserve fMRI signal details, and eventually create a general CS framework.

In this study, we tested the feasibility of CS for non-EPI fMRI at 9.4T using the model of rat somatosensory stimulation. Fully sampled data with high-resolution spoiled gradient echo and 4 independent pass-band bSSFP fMRI each with different phase-cycling angles (0°, 90°, 180°, and 270°) underwent retrospective downsampling and reconstruction using k-t FOCUSS algorithm. Various sampling patterns and sparsifying transforms such as temporal FT and KLT were employed to systematically study the effects of different k-space sampling pattern and the effects of choosing a different CS reconstruction algorithm in high field CS fMRI. The baseline image quality and sensitivity and specificity of activation maps from data with CS reconstruction were compared to those from the original full-sampled data. The potential for improving the temporal resolution of non-EPI fMRI at high magnetic fields without sacrificing quality of fMRI activation maps is demonstrated in this paper.

2. Methods

2.1. Animal Preparation and Data Acquisition

Three male Sprague-Dawley rats weighing 250~450 g (Charles River Laboratories, Wilmington, MA, USA) were used with the approval from the Institutional Animal Care and Use Committee (IACUC) at University of Pittsburgh. Animal preparation was the same as previously published [7]. Briefly, the rats were intubated for mechanical ventilation (RSP-1002, Kent Scientific, CT, USA). The catheters were inserted in the femoral artery and femoral vein for blood gas sampling and fluid administration (5% dextrose in saline infused at 0.4 mL/hr), respectively. Once the surgery was finished, the isoflurane level was maintained at 1.4%. Ventilation rate and volume were adjusted based on blood gas analysis results (Stat profile pHOx; Nova Biomedical, MA, USA).

Electrical stimulation was applied to either the right or left forepaw using two needle electrodes inserted under the skin between digits 2 and 4 [23]. Stimulation parameters for activation studies were as follows: current = 1.2~1.6 mA, pulse duration = 3 ms, repetition rate = 6 Hz, stimulation duration = 15 s, and interstimulation period = 3 min.

All experiments were carried out on a Varian 9.4T/31 cm MRI system (Palo Alto, CA) with an actively shielded gradient coil of 12 cm inner diameter, which operates at a maximum gradient strength of 40 G/cm and rise time of 120 μs. A homogeneous coil and a surface coil (Nova Medical, Wilmington, MA) were used for RF excitation and reception, respectively. Localized shimming was performed with point resolved spectroscopy [24] over a coronal slab (~12 × 6 × 6 mm3) covering forelimb somatosensory cortex to yield a water spectral linewidth of 20~30 Hz. Spoiled gradient echo (which is denoted by GRE throughout this paper) and pass-band bSSFP studies were performed with TR/TE = 20/10 ms and 10/5 ms, respectively. The bSSFP fMRI studies were performed with four different phase-cycling angles (θ) of 0°, 90°, 180°, and 270° (which are denoted by PC 0, PC 1, PC 2, and PC 3, resp., throughout this paper). The resolution parameters were the same for all studies: matrix size = 256 × 192, FOV = 2.4 × 2.4 cm2, number of slice = 1, and slice thickness = 2 mm. Flip angles for all the bSSFP and GRE fMRI studies were 16° and 8°, respectively. Forty-eight measurements were acquired for each bSSFP fMRI study: 16 during prestimulus baseline, 8 during stimulation, and 24 during the poststimulus period. These numbers of measurements were reduced by half for GRE fMRI study, in order to maintain the same spatial resolution. Four bSSFP and one GRE fMRI studies composed one full set and each full set was repeated 15 to 25 times for averaging per subject rat.

2.2. k-Space Sampling Patterns

The initial data to undergo CS reconstruction is important to guarantee high performance of CS algorithms. All fMRI studies in this work were conducted with block design paradigm; thus downsampling was considered in the k-t (i.e., k-space-temporal) domain to utilize CS algorithms optimized for exploiting temporal redundancy in dynamic MRI. In order to determine the optimal sampling pattern for CS application on fMRI data acquired at high field, four different sampling patterns were considered (Figures 1(b)–1(e)): sampling masks were generated using uniform random sampling (Figure 1(b)), Gaussian random sampling (Figure 1(c)), a mixture of Gaussian and uniform random sampling (Figure 1(d)), and a mixture of Gaussian and uniform random sampling with full sampling of k-space center 1 line (Figure 1(e)). The generated sampling masks were applied to each full-sampled k-space dataset for retrospective downsampling before the CS reconstruction procedure.

Figure 1

Masks of five different sampling patterns with downsampling factor of 4. Patterns of (a) full sampling of k-space center, (b) uniform random sampling, (c) Gaussian random sampling, (d) mixture of Gaussian and uniform random sampling, and (e) mixture of Gaussian and uniform random sampling with full sampling of center 1 line are shown. All masks are generated to sample only a quarter of the total data (total size of the data indicated by dashed box). Notice that no k-space high frequency information is included in (a) and (c).

In this study, a fixed downsampling factor of 4 was applied to all datasets: only a quarter of the original k-space data was used for CS reconstruction. Each 2D sampling mask was generated for each time frame and a fixed number of N PE/4 lines were sampled along the phase-encoding (PE) direction to maintain the downsampling factor of 4 (where N PE indicates the total number of PE lines). Note that bSSFP or GRE sequences utilize different RF excitations for each TR; thus the acceleration factor depends on the total number of PE lines sampled. The uniform random sampling mask was generated by sampling the PE lines according to a uniform probability distribution. The Gaussian random sampling mask was generated by sampling the PE lines according to a Gaussian probability distribution of P(k ) = Ae −(, with σ = N PE/9 (where k indicates k-space PE line number index in the range of −95 to 96, σ indicates standard deviation, and A indicates the weighting factor which was adjusted to make the equation a valid probability density function). The mixture of Gaussian and uniform random sampling mask was generated by sampling N PE/6 number of PE lines according to the above Gaussian probability distribution and subsequently sampling N PE/12 number of PE lines according to a uniform random probability distribution (N PE/6 + N PE/12 = N PE/4 lines were sampled in total). The mixture of Gaussian and uniform random sampling with full sampling of k-space center 1 line was generated similarly, but with continuous sampling of the k-space center 1 line (i.e., k = 0) for each time frame along the temporal dimension. Pseudocodes for generation of the sampling patterns are provided as follows.

Pseudocode for Generation of k-Space Downsampling Pattern on 2D Time-Series MRI Data.

For downsampling factor (D), one has the following.

Uniform Random Sampling Mask

Generate uniform probability distribution along the PE dimension of k-space data.

Sample N PE/D number of PE lines according to the uniform probability distribution.

Repeat steps (1)-(2) for each time frame.

Gaussian Random Sampling Mask

Define σ in relation to N PE.

Generate Gaussian probability distribution P(k ) = Ae −( along the PE dimension of k-space data.

Sample N PE/D number of PE lines according to the Gaussian probability distribution.

Repeat steps (1)–(3) for each time frame.

Mixture of Gaussian and Uniform Random Sampling Mask (Gaussian : Random = a : b)

Define σ in relation to N PE.

Generate Gaussian probability distribution P(k ) = Ae −( along the PE dimension of k-space data.

Sample (N PE/D)×(a/(a + b)) number of PE lines according to the Gaussian probability distribution.

Generate uniform probability distribution along the PE dimension for the remaining PE lines.

Sample (N PE/D)×(b/(a + b)) number of PE lines according to the uniform probability distribution.

Repeat steps (1)–(5) for each time frame.

Mixture of Gaussian and Uniform Random Sampling Mask with Full Sampling of Center 1 Line (Gaussian : Random = a : b)

Sample 1 PE line from the k-space center (e.g., k = 0).

Repeat all steps from (C) with N PE/D − 1 number of PE lines.

The Gaussian probability distribution was considered as a derivative form of k-space center-weighted downsampling pattern, with stronger weighting on k-space low-frequency information. Sampling using a mixture of Gaussian and uniform random probability distributions was considered as a further variation preserving information at both high and low frequencies. The inclusion of k-space center 1 line was considered as an option to preserve the lowest frequency information, since the center of k-space contains most of the contrast information for an MR image. For comparison analysis, original full-sampled data was used as the ground truth to study the effect of CS reconstruction, and a simple quarter downsampled mask with full sampling of k-space center (Figure 1(a)) was used as a control. Throughout the paper, uniform random sampling pattern, Gaussian random sampling pattern, mixture of Gaussian and uniform random sampling pattern, and mixture of Gaussian and uniform random sampling pattern with full sampling of center 1 line will be denoted by PatternR, PatternG, PatternGR, and PatternGRC1, respectively.

2.3. CS Algorithms

Two variations of k-t FOCUSS algorithms, temporal FT and KLT (which will be further denoted by Algorithms 1 and 2, resp.), were used in this study (see the Appendix for detailed description of k-t FOCUSS algorithms) [16]. The covariance matrix for Algorithm 2 was defined to be constructed from an initial reconstruction using Algorithm 1 with preliminary parameter of N FOC1 = 2 for each dataset:where indicates the reconstruction from Algorithm 1 (please refer to the Appendix for detailed definition of ). The eigenvectors of the covariance matrix were further used as the KL transform (Φ) and N KLT = 1 was used to update Φ once.

Algorithm 1

= ktFOCUSS(p 1, λ 1, N CG1, N FOC1, Φ, x (0)).

Algorithm 2

k-t FOCUSS with Karhunen-Loéve transform.

2.4. k-t FOCUSS Parameters

For both algorithms, weighting matrix power factor (p) of 0.5 was used to find the sparse solution equivalent to the l 1 solution of CS [16]. Conjugate gradient (CG) iteration number (N CG) of 30 was considered sufficient and was used for both Algorithm 1 (N CG1) and Algorithm 2 (N CG2), based on previous application with k-t FOCUSS [25]. Regularization factor (λ) of 0.1 was used for Algorithm 1 (λ 1) [16] and 0.01 was used for Algorithm 2 (λ 2) [22].

Both Algorithms 1 and 2 were optimized based on variation in the FOCUSS iteration number (N FOC) parameters (i.e., N FOC1 and N FOC2, resp.). The following stopping criterion is used to determine the optimal N FOC value from each dataset in training phase:where denotes the CS reconstruction of the spatiotemporal fMRI data at kth FOC iteration, denotes the CS reconstruction of the spatiotemporal fMRI data at (k − 1)th FOC iteration, and ‖·‖ denotes the Frobenius norm. The performance of k-t FOCUSS algorithms with the proposed stopping criterion for N FOC parameters was evaluated via subject-based leave-one-out cross-validation (i.e., the N FOC value for data from each subject was determined based on a training set consisting of data from the other subjects).

The effects of N CG and λ were investigated separately for verification with the determined optimal N FOC value (i.e., value found during the training phase) for each subject data. The reconstruction of phase-cycled bSSFP and GRE data was used for investigation with all cases of sampling patterns as follows. Residual error = was used as a measure to observe data fitting and convergence in each k-t FOCUSS algorithm with increase in N CG, where y denotes the k-space time-series data of the sampled k-space lines and denotes the CS reconstruction of the k-space time-series data for the corresponding k-space lines. Average mean square error (MSE) of the whole time-series data was used as a measure to observe the noise level in the reconstructed image with variation in λ, which was calculated as follows:where u denotes the original full-sampled spatiotemporal fMRI data at time frame t, denotes the CS reconstructed spatiotemporal fMRI data at time frame t, T denotes the total number of time frames, and N denotes the total number of image pixels at each time frame.

2.5. Region of Interest Selection

A region of interest (ROI) was selected to help compare the effects of different sampling patterns and CS algorithms. The regions determined to be functionally active (i.e., rejecting the null hypothesis H 0) according to the t-statistics map of the original full-sampled data were chosen as the ROI for further analysis. New ROIs were defined for each dataset.

2.6. Quantitative Analysis

Frame-by-frame normalized MSE, t-statistics functional map, ROI time course plot, and receiver operating characteristics (ROC) curve were calculated for further investigation of the applicability of CS for fMRI data at high field. These analyses were performed on the bSSFP data with PC2 for clear evaluation of the effect of CS application, since the sequence corresponds to the conventional bSSFP sequence (i.e., 180° phase-cycling) displaying a fairly uniform signal contrast without any significant banding artifacts and showed clear activation foci in the full-sampled data.

The frame-by-frame normalized MSE calculation at time t was performed using the following equation:Student's t-test was performed for each dataset to statistically analyze fMRI data and generate the t-statistics functional map. The T-score is calculated on a pixel by pixel basis over time as follows:where denote the mean, s · denote the standard deviation, and n · denote the length of the baseline time-series x and activation time-series y, respectively. The t-statistics functional map was generated for a significance level of 0.05, and clusters less than 6 pixels were rejected. ROI time course was plotted as the mean ROI value.

The ROC curve was generated to provide standardized and statistically meaningful means for comparing fMRI signal-detection accuracy [26]. For each dataset, the t-statistics map generated from the original full-sampled data with significance level of 0.05 was used as the ground truth. True positive fraction (TPF) and false positive fraction (FPF) were calculated over various significance levels to generate the ROC curve. The performance was measured by the area under the curve (AUC) ranging from 0 to 1, with 1 representing better performance. The TPF and FPF were calculated using the following equations:where TPF relates to sensitivity and 1 − FPF relates to specificity.

3. Results

3.1. Determination of N FOC

N FOC values determined via subject-based leave-one-out cross-validation for Algorithm 1 were 4, 3, 3, and 3 for PatternR, PatternG, PatternGR, and PatternGRC1, respectively, and those for Algorithm 2 were 5, 4, 4, and 4 for PatternR, PatternG, PatternGR, and PatternGRC1, respectively. Identical N FOC values were found regardless of the pulse sequence type (i.e., GRE and bSSFP PC0, PC1, PC2, and PC3) in all the subjects. All further analyses were performed with the determined N FOC values for each subject data, to evaluate the performance of k-t FOCUSS algorithms with the proposed stopping criterion.

3.2. Original Data: High Field bSSFP

Different phase-cycling angles in the fMRI maps of full-sampled bSSFP data showed shifting in activation foci (i.e., the activation foci were located around the cortical surface area for PC1 and 2, while they were located in the middle cortical regions for PC0 and 3, as indicated by white arrows in Figure 2(b)). This spatial shift of activation foci as a function of PC angle implies that the high field phase-cycled bSSFP maps are spatially heterogeneous due to magnetic field inhomogeneity and was used in this study to confirm that CS with an appropriate downsampling scheme can preserve the details of the spatial pattern of the functional activation.

Figure 2

Baseline and fMRI images of original and control. (a) Baseline images and (b) fMRI maps of full-sampled data and (c) baseline images and (d) fMRI maps of downsampled data with full acquisition of k-space center are shown. The 20th and 10th time frame of bSSFP and GRE is shown for baseline images, respectively. The fMRI maps are shown for significance level of 0.05. Downsampling pattern and acquisition method are shown on the top and left-hand side of the images, respectively. The DC artifact is indicated by yellow arrow and the enlarged fMRI activation region is indicated by blue dashed box. Notice the activation foci shift (white arrow) in the phase-cycled bSSFP data. The case from a representative rat is shown since similar results were obtained for different subject rats.

3.3. Reconstruction with Algorithm 1: k-t FOCUSS with Temporal FT

Visually the original baseline images became blurred with artifacts after downsampling was applied (Figure 3). Despite the distortion and degradation after downsampling, the baseline images were well reconstructed using Algorithm 1 regardless of sampling pattern (Figures 4(a)–4(d)). Visually the image contrast and resolution were well preserved for CS reconstructed images from all sampling patterns compared to the original baseline image (Figure 2(a)) and downsampled baseline image with only k-space low-frequency information (Figure 2(c)). The frame-by-frame normalized MSE values from all the CS sampling patterns were significantly lower than those from downsampling with only k-space low-frequency information (Figure 5), indicating high reconstruction performance of Algorithm 1. In particular, the mixture of Gaussian and uniform random sampling scheme (i.e., PatternGR and PatternGRC1) showed the lowest frame-by-frame normalized MSE values across all time frames. Overall, all Gaussian-weighted sampling patterns showed increased spatial resolution and SNR (Figures 4(b), 4(c), and 4(d)) with reduction of artifacts (indicated by yellow arrow in Figure 2(a)) for all phase-cycled bSSFP data, while downsampling with only k-space low-frequency information showed increase of artifacts (Figure 2(c)).

Figure 3

Baseline images after downsampling. Baseline images of downsampled data using (a) PatternR, (b) PatternG, (c) PatternGR, and (d) PatternGRC1 are shown. The 20th and 10th time frame of bSSFP and GRE is shown, respectively. Downsampling pattern and acquisition method are shown on the top and left-hand side of the images, respectively. Notice the blurring and artifacts in the downsampled images. The case from a representative rat is shown since similar results were obtained for different subject rats.

Figure 4

Comparison of baseline images reconstructed using Algorithm 1 (k-t FOCUSS with temporal FT). Baseline images of CS reconstructed data using Algorithm 1 and (a) PatternR, (b) PatternG, (c) PatternGR, and (d) PatternGRC1 are shown. The 20th and 10th time frame of bSSFP and GRE is shown, respectively. Downsampling pattern and acquisition method are shown on the top and left-hand side of the images, respectively. The case from a representative rat is shown since similar results were obtained for different subject rats.

Figure 5

Comparison of frame-by-frame normalized MSE plots from reconstructed baseline images using Algorithm 1 (k-t FOCUSS with temporal FT). The frame-by-frame normalized MSE plots of bSSFP time-series data with PC angle of 180° are shown. The case from a representative rat is shown since similar results were obtained for different subject rats.

The fMRI maps were also reconstructed well from all Gaussian-weighted sampling patterns using Algorithm 1 (Figures 6(b), 6(c), and 6(d)), while those from PatternR did not show any meaningful functional activations (Figure 6(a)). The fMRI maps from downsampling with only k-space low-frequency information showed significant blurring in the activation region (Figure 2(d)). The fMRI maps from PatternGR (Figure 6(c)) and PatternGRC1 (Figure 6(d)) were closer to the original fMRI maps than those from PatternG (Figure 6(b)) in terms of preserving details in activation foci shift, presumably due to the inclusion of appropriate high frequency k-space information.

Figure 6

Comparison of fMRI maps reconstructed using Algorithm 1 (k-t FOCUSS with temporal FT). The fMRI maps of CS reconstructed data using Algorithm 1 and (a) PatternR, (b) PatternG, (c) PatternGR, and (d) PatternGRC1 are shown for significance level of 0.05. Downsampling pattern and acquisition method are shown on the top and left-hand side of the images, respectively. Notice the activation foci shift (white arrow) in the phase-cycled bSSFP data. The case from a representative rat is shown since similar results were obtained for different subject rats.

The time course of the mean ROI value was also relatively well preserved in images from all Gaussian-weighted sampling patterns (Figures 7(b), 7(c), and 7(d)), while those from PatternR differed from the original with significantly increased temporal fluctuation (Figure 7(a)). Mean ROI time courses of images from all Gaussian-weighted sampling patterns resembled those of the original data; even with only (1/4)th of the whole data, the mean ROI time course followed the trend of the original time course with slightly reduced mean amplitude difference, percent signal change, and also signal fluctuation. These observations were applicable regardless of acquisition method and different PC angles for bSSFP. The AUC value of ROC curves indicated overall high sensitivity and specificity of all Gaussian-weighted sampling patterns using Algorithm 1 (Table 1). PatternGRC1 displayed the highest ROC performance than other sampling patterns including downsampling with only k-space low-frequency information.

Figure 7

Comparison of mean ROI time course plots between full-sampled original data and CS reconstructed data using Algorithm 1 (k-t FOCUSS with temporal FT). The time course of mean ROI from CS reconstructed data using Algorithm 1 and (a) PatternR, (b) PatternG, (c) PatternGR, and (d) PatternGRC1 is shown. Mean ROI is calculated from bSSFP time-series data with PC angle of 180°. The case from a representative rat is shown since similar results were obtained for different subject rats.

Table 1

ROC curve AUC values from fMRI maps of bSSFP time-series data with PC angle of 180° of a representative rat.

4x Down Center (control)		CS
		PatternR	PatternG	PatternGR	PatternGRC1
0.9794	Algorithm 1	0.8943	0.9800	0.9795	0.9827
	Algorithm 2	0.8639	0.9726	0.9714	0.9825

The reconstruction times of the k-t FOCUSS algorithms are shown for bSSFP PC2 and GRE in Table 2. Only one representative case is shown for each bSSFP and GRE data since similar results were obtained regardless of bSSFP PC angle, sampling pattern, and subject rats, despite the difference in total time due to the usage of different N FOC parameter values. The reconstruction time for the bSSFP data was approximately twice as long as that of the GRE data, since the speed of reconstruction is largely dependent on the size of the matrix and number of time frames (recall that the number of time frames of GRE data was half of that of bSSFP data).

Table 2

Reconstruction speed of k-t FOCUSS algorithms on downsampled time-series data using PatternGRC1 of a representative rat.

	Algorithm 1		Algorithm 2
	Single FT-CG iteration	Total time	Covariance calculation	Single KLT-CG iteration	Total time
bSSFP PC2	0.46 s	41.35 s	27.81 s	0.64 s	77.26 s
GRE	0.24 s	21.29 s	14.27 s	0.29 s	34.98 s

FT: Fourier transform, CG: conjugate gradient method, and KLT: Karhunen-Loéve transform.

Calculated on PC (Windows 7), CPU: 3.30 GHz, RAM: 4.00 GB.

3.4. Reconstruction with Algorithm 2: k-t FOCUSS with KLT

Results of Algorithm 2 in terms of baseline images, frame-by-frame normalized MSE plot, fMRI maps, ROI time course, and AUC values in ROC curve are shown in Figures 8, 9, 10, and 11 and Table 1, respectively. Overall the results were similar to those of Algorithm 1. Slight differences were observed between algorithms in the view points of sensitivity and specificity depending on sampling pattern. The images from all sampling patterns for Algorithm 2 showed slightly lower sensitivity and specificity than those for Algorithm 1. The reconstruction time of Algorithm 2 was longer than that of Algorithm 1, mainly due to the calculation of covariance matrix requiring the preliminary estimation (Table 2).

Figure 8

Comparison of baseline images reconstructed using Algorithm 2 (k-t FOCUSS with KLT). Baseline images of CS reconstructed data using Algorithm 2 and (a) PatternR, (b) PatternG, (c) PatternGR, and (d) PatternGRC1 are shown. The 20th and 10th time frame of bSSFP and GRE is shown, respectively. Downsampling pattern and acquisition method are shown on the top and left-hand side of the images, respectively. The case from a representative rat is shown since similar results were obtained for different subject rats.

Figure 9

Comparison of frame-by-frame normalized MSE plots from reconstructed baseline images using Algorithm 2 (k-t FOCUSS with KLT). The frame-by-frame normalized MSE plots of bSSFP time-series data with PC angle of 180° are shown. The case from a representative rat is shown since similar results were obtained for different subject rats.

Figure 10

Comparison of fMRI maps reconstructed using Algorithm 2 (k-t FOCUSS with KLT). fMRI maps of CS reconstructed data using Algorithm 2 and (a) PatternR, (b) PatternG, (c) PatternGR, and (d) PatternGRC1 are shown for significance level of 0.05. Downsampling pattern and acquisition method used are shown on the top and left-hand side of the images, respectively. Notice the activation foci shift (white arrow) in the phase-cycled bSSFP data. The case from a representative rat is shown since similar results were obtained for different subject rats.

Figure 11

Comparison of mean ROI time course plots between full-sampled original data and CS reconstructed data using Algorithm 2 (k-t FOCUSS with KLT). The time course of mean ROI from CS reconstructed data using Algorithm 2 and (a) PatternR, (b) PatternG, (c) PatternGR, and (d) PatternGRC1 is shown. Mean ROI is calculated from bSSFP time-series data with PC angle of 180°. The case from a representative rat is shown since similar results were obtained for different subject rats.

3.5. Investigation of N CG and Effect

Representative results from CS reconstruction of PatternGRC1 using Algorithms 1 and 2 are shown in Figures 12, 13, 14, and 15, and similar results were obtained regardless of acquisition method and sampling patterns. Based on the results from Figures 12 and 14, the N CG value of 30 (i.e., value used for both Algorithms 1 and 2 throughout the paper) seems to be sufficient to ensure data fitting and error convergence. As shown in Figures 13 and 15, the reconstruction error of k-t FOCUSS algorithms was relatively insensitive to λ variation and minimal error was achieved with small values of λ (e.g., less than 1).

Figure 12

Example showing the effect of N CG1 on reconstruction using Algorithm 1 (k-t FOCUSS with temporal FT). Residual error with N CG1 variation using Algorithm 1 on (a) bSSFP PC0, (b) bSSFP PC1, (c) bSSFP PC2, (d) bSSFP PC3, and (e) GRE data is shown for each different N FOC1. Plots are obtained from reconstruction of downsampled data using PatternGRC1 during testing phase of a representative rat, and similar results were obtained from different sampling patterns and different subject rats. Notice that residual error reaches convergence as N CG1 increases. The chosen N CG1 value of 30 is represented in each plot as a red dashed line.

Figure 13

Example showing the effect of λ 1 on reconstruction using Algorithm 1 (k-t FOCUSS with temporal FT). Average MSE with λ 1 variation using Algorithm 1 on (a) bSSFP PC0, (b) bSSFP PC1, (c) bSSFP PC2, (d) bSSFP PC3, and (e) GRE data is shown. Plots are obtained from reconstruction of downsampled data using PatternGRC1 and optimal N FOC1 found during testing phase of a representative rat, and similar results were obtained from different sampling patterns and different subject rats. Notice that reconstruction error is relatively invariant to λ 1 variation and minimal error is achieved with small values of λ 1.

Figure 14

Example showing the effect of N CG2 on reconstruction using Algorithm 2 (k-t FOCUSS with KLT). Residual error with N CG2 variation using Algorithm 2 on (a) bSSFP PC0, (b) bSSFP PC1, (c) bSSFP PC2, (d) bSSFP PC3, and (e) GRE data is shown for each different N FOC2. Plots are obtained from reconstruction of downsampled data using PatternGRC1 during testing phase of a representative rat, and similar results were obtained from different sampling patterns and different subject rats. Notice that residual error reaches convergence as N CG2 increases. The chosen N CG2 value of 30 is represented in each plot as a red dashed line.

Figure 15

Example showing the effect of λ 2 on reconstruction using Algorithm 2 (k-t FOCUSS with KLT). Average MSE with λ 2 variation using Algorithm 2 on (a) bSSFP PC0, (b) bSSFP PC1, (c) bSSFP PC2, (d) bSSFP PC3, and (e) GRE data is shown. Plots are obtained from reconstruction of downsampled data using PatternGRC1 and optimal N FOC2 found during testing phase of a representative rat, and similar results were obtained from different sampling patterns and different subject rats. Notice that reconstruction error is relatively invariant to λ 2 variation and minimal error is achieved with small values of λ 2.

4. Discussion

To our knowledge, it is the first study to apply CS to high field bSSFP fMRI and to systematically evaluate effects of CS sparsity schemes on non-EPI fMRI. The CS sampling scheme should be determined in relation to the CS algorithm to preserve detailed image information appropriately. Dense sampling in k-space center region such as Gaussian-weighting or inclusion of k-space center 1 line was better for CS, due to the fact that most energy is located in the k-space center region and also due to incoherent aliasing effects from variable-density sampling [13, 16]. The reconstruction results from PatternGR and PatternGRC1 also verify that a variation in sampling scheme with more k-space high frequency information leads to better reconstruction performance preserving signal details (e.g., activation foci shifting phenomenon in bSSFP), which may become critical for applications such as fMRI studies. Inclusion of more k-space low-frequency information implies less k-space high frequency information which may lead to an enlargement or blurring of the activation foci in fMRI maps (e.g., Figures 2(d), 6(b), and 10(b)). Thus, both k-space center and edge regions are important, and methods that achieve a certain balance between them need to be exploited for correct reconstruction of non-EPI fMRI data using CS. Overall the mixture of Gaussian and uniform random sampling scheme reconstructed both the baseline images and fMRI maps well while preserving the signal details and thus seems to be an ideal sampling scheme for CS applied to non-EPI fMRI.

The two algorithms of k-t FOCUSS with temporal FT and KLT showed similar performances overall. The slight differences in their results are presumed to be due to the utilization of different transformation domains for each iteration. Interestingly, k-t FOCUSS with temporal FT performed slightly better than k-t FOCUSS with KLT in terms of ROC performance in this study, despite the fact that KLT is known as an efficient spectral decorrelator [27, 28]. Several factors may account for this. First, the fMRI studies were performed with block design paradigm in this work, and temporal redundancy from the spatial-temporal frequency domain may have been exploited better for such data type. Since KLT is a data-driven transform, k-t FOCUSS with KLT may potentially perform better in cases of rapid event-related paradigms. Second, the decorrelation of nonperiodic noise might not have been noticeable in the image, since bSSFP sequence is known to provide the highest SNR per unit time [29, 30] and the simulation studies were performed on datasets with enough averaging (e.g., 15 to 25 times). Recently, it has been reported that application of CS to fMRI can increase FPF in real acquisition settings, and k-t FOCUSS with KLT has shown to reconstruct fMRI maps with reduced false activations [22]. Thus, the effectiveness of both algorithms needs verification with real fMRI studies. Nonetheless, results from the current study indicate that both algorithms are potentially good solutions for acceleration of high field non-EPI fMRI.

Appropriate choice of CS reconstruction parameters is one of the main concerns of the application of CS. The optimal parameters may vary depending on noise level, temporal resolution, and other possible factors in actual data acquisition environment. In general, reconstruction parameters are found with known noise level [31] or alternatively are selected via cross-validation [32-34]. There are multiple parameters involved for the case of k-t FOCUSS algorithm, which requires hyperparameter optimization and thus increases the computation burden [16]. Therefore, the effect of two different k-t FOCUSS parameters, λ and N CG, is additionally investigated in this study. Considering the physical meaning of each parameter, two different metrics were used for evaluation. Since the regularization parameter λ is a tuning parameter used to find the solution with best improvement in SNR, average MSE was used to show its effect on the noise level in the reconstructed image. Since the CG method is employed to iteratively find the solution to the unknown signal (i.e., denoted by x in (A.9) of the Appendix), residual error was used to investigate data fitting and convergence with decreases in measurement error (i.e., difference between sampled k-space measurements y and estimation ) as number of CG iterations increases (note that the signal-measurement relationship is defined in (A.4) of the Appendix and can be used to find from x). Based on the results from Figures 12 and 14, fixation of N CG to a sufficient value (e.g., 30 in case of our study) and λ to a small value is preferred to reduce parameter variability and to simplify the usage of k-t FOCUSS algorithms on high field non-EPI fMRI studies. These results agree with previous applications of k-t FOCUSS where a sufficient value of N CG and a small value of λ are used and are proven to perform well in high-quality fMRI studies from real scanner acquisitions [16, 22]. However, care should be taken for extrapolation of these parameters for data types different from those of the current studies. With fixation of N CG and λ, the only issue for the application of k-t FOCUSS algorithms for high field bSSFP fMRI data lies in the choice of N FOC. The N FOC parameters found from the current study need to be tested in the context of real CS application for verification and general usage. The choice of a high N FOC value may ensure minimal error for most cases of applications; however, this also leads to increased number of calculations required for CS reconstruction. Thus, the trade-off between minimization of error and increase in postprocessing time must be considered appropriately before choosing the N FOC value for future studies.

Eddy currents can cause problems in bSSFP imaging with nonlinear phase-encoding orders. The problems may become noticeable when the sparsity schemes tested in this study are implemented in real acquisition settings. Previously Bieri et al. [35] discovered a simple method to suppress the eddy current effect by pairing two consecutive k-space lines. By incorporating this idea, a simulation study was performed to see the effect of pairwise downsampling scheme. A paired PatternGRC1 was generated with a downsampling factor of 4. The downsampling scheme and reconstructed fMRI maps from each k-t FOCUSS algorithm are shown in Figure 16, and the ROC performance of the reconstructed data is shown in Table 3. The employment of pairwise sampling scheme showed maintenance of activation foci shift but decrease in both activation detection sensitivity and specificity compared to the results without pairwise sampling (Figures 6(d) and 10(d)), regardless of PC angle and k-t FOCUSS algorithm. Overall, these results indicate that the pairwise sampling scheme may be used to suppress eddy current artifacts in bSSFP fMRI with CS, but there exists trade-off between the suppression of eddy current effect and fMRI sensitivity as well as specificity.

Figure 16

Presumed effect of paired random sampling scheme in the reconstruction of fMRI map using k-t FOCUSS algorithms. The fMRI maps of CS reconstructed data using (a) Algorithm 1 and pairwise sampling of PatternGRC1 and (b) Algorithm 2 and pairwise sampling of PatternGRC1 are shown for significance level of 0.05. Downsampling pattern and acquisition method used are shown on the top and left-hand side of the images, respectively. Notice the decrease in activation sensitivity with the inclusion of pairwise sampling scheme.

Table 3

ROC curve AUC values from reconstructed fMRI maps using k-t FOCUSS algorithms and pairwise sampling of PE lines with PatternGRC1 on bSSFP time-series data with PC angle of 180° of a representative rat.

	Algorithm 1	Algorithm 2
Pairwise PE PatternGRC1	0.9810	0.9736

Application of CS to high field non-EPI fMRI can be meaningful for high-resolution fMRI studies, since conventional GRE-EPI fMRI is sensitive to image distortion and degradation caused by local magnetic field inhomogeneity at high magnetic fields. Although the temporal resolution of non-EPI sequences is lower than the typically used GRE-EPI, it is shown through the study that the temporal resolution or the spatial coverage can be improved using CS. Several potential advantages of CS can be derived for fMRI studies in this regard. First, better temporal resolution increases the number of time frames within a given time and can in turn improve the statistical power of BOLD activations [22, 36]. Second, the weighted-norm process of the CS algorithm can reduce artifacts from scanner-related drifts, respiratory-induced noise, cardiac pulsation, and subject motion [37-39] and can also improve the activation detection sensitivity in t-statistics [22]. Lastly, CS can improve spatial coverage which is essential for many fMRI studies that require a big ROI or ROIs from multiple brain regions. Thus, the application of CS in fMRI has great potential in practice.

One negative aspect of CS in fMRI studies is the addition of postprocessing time related to CS reconstruction. The reconstruction time of k-t FOCUSS algorithms is largely affected by the iteration parameters and the size of the data (e.g., matrix size, number of slices, number of time frames, etc.). The results from Table 2 imply that the temporal resolution or spatial coverage of non-EPI sequence fMRI studies can be improved using CS at the cost of reasonable addition of postprocessing time (i.e., several minutes). Therefore, depending on applications, the trade-off between reconstruction time and temporal resolution (and/or spatial coverage) must be investigated before applying CS algorithms to fMRI studies.

In the past, there have been many dynamic MRI studies other than fMRI with acceleration factor of 8 or higher using CS [18, 40–42]. However, CS has been applied to fMRI in a limited number of studies and acceleration factor up to 4 was used in most of the truly accelerated fMRI studies [21, 22]. Decrease in image quality has also been reported in some pilot studies after CS reconstruction even with 2-fold acceleration [21]. This may be attributed to the fact that distinct from other dynamic MRI studies fMRI requires detection of fine signal changes, which can be achieved by preserving high frequency information. Based on the results from our studies, acceleration factor of 4 seems sufficient for CS application on high field non-EPI fMRI studies. For example, for a bSSFP fMRI experiment with matrix size = 128 × 128 and TR = 5 ms, the temporal resolution becomes 0.64 s for a single slice. Thus, the 4-fold acceleration can improve the temporal resolution up to 0.16 s (i.e., close to the temporal resolution of EPI) or the spatial coverage up to 24 slices (i.e., near whole brain coverage) with temporal resolution less than 4 s (note that although TR was 10 ms in this study, bSSFP with TR ≤ 10 ms has been successfully applied to fMRI at high field ≥7T). Nonetheless, improvements can be made to increase the acceleration factor above 4 and the quality may depend on the scan condition (e.g., scan resolution, SNR). A systematic study for different acceleration factors may need to be conducted to further improve application of CS on high field non-EPI fMRI studies.

In this paper, the effect of CS is investigated through retrospective downsampling of full-sampled fMRI data. Therefore, further works are necessary to verify the downsampling schemes that were retrospectively evaluated in this study, by implementing them and performing real fMRI studies, which is beyond the scope of this paper.

5. Conclusion

The CS reconstruction of fMRI data acquired at high field using k-t FOCUSS varies greatly with sampling scheme and thus the sampling scheme must be selected appropriately. Information in both k-space low and high frequency regions is important for better reconstruction performance and preservation of signal details, respectively, and thus sampling schemes that achieve a certain balance between the two must be selected for the application of CS to non-EPI fMRI data. The two k-t FOCUSS algorithms, temporal FT and KLT, showed good reconstruction results overall with effective suppression of downsampling artifacts and improved spatial resolution and thus are good candidates for CS in high field non-EPI fMRI studies. The application of CS to fMRI has great potential in practice for improvement of temporal resolution and/or spatial coverage.