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.4T 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.


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][3][4][5][6][7][8][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][11][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 -space filling trajectory. The other solution to improve the temporal resolution of non-EPI sequences is to use compressed sensing (CS) [13][14][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 -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 -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.

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 mm 3 ) 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 cm 2 , number of slice = 1, and slice thickness = 2 mm. Flip angles for all the bSSFP and GRE fMRI studies were 16 ∘ and 8 ∘ , respectively. Fortyeight 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.

-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 -(i.e., -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 -space center 1 line (Figure 1(e)). The generated sampling masks were applied to each full-sampled -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 -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 -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 -space data was used for CS reconstruction. Each 2D sampling mask was generated for each time frame and a fixed number of PE /4 lines were sampled along the phase-encoding (PE) direction to maintain the downsampling factor of 4 (where 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 ( ) = −( ) 2 /2 2 , with = PE /9 (where indicates -space PE line number index in the range of −95 to 96, indicates standard deviation, and 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 PE /6 number of PE lines according to the above Gaussian probability distribution and subsequently sampling PE /12 number of PE lines according to a uniform random probability distribution ( PE /6 + PE /12 = PE /4 lines were sampled in total). The mixture of Gaussian and uniform random sampling with full sampling of -space center 1 line was generated similarly, but with continuous sampling of the -space center 1 line (i.e., = 0) for each time frame along the temporal dimension. Pseudocodes for generation of the sampling patterns are provided as follows.

Pseudocode for Generation of -Space Downsampling Pattern on 2D Time-Series MRI Data.
For downsampling factor ( ), one has the following. (1) Define in relation to PE .
(2) Repeat all steps from (C) with PE / −1 number of PE lines.
The Gaussian probability distribution was considered as a derivative form of -space center-weighted downsampling pattern, with stronger weighting on -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 -space center 1 line was considered as an option to preserve the lowest frequency information, since the center of -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 -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 Pattern R , Pattern G , Pattern GR , and Pattern GRC1 , respectively.

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 FOC1 = 2 for each dataset: whereÛ = [û 1 ,û 2 , . . . ,û ] 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 KLT = 1 was used to update Φ once.

k-t FOCUSS Parameters.
For both algorithms, weighting matrix power factor ( ) of 0.5 was used to find the sparse solution equivalent to the 1 solution of CS [16]. Conjugate gradient (CG) iteration number ( CG ) of 30 was considered sufficient and was used for both Algorithm 1 ( CG1 ) and Algorithm 2 ( 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 ( FOC ) parameters (i.e., FOC1 and FOC2 , resp.). The following stopping criterion is used to determine the optimal FOC value from each dataset in training phase: whereÛ ( ) denotes the CS reconstruction of the spatiotemporal fMRI data at th FOC iteration,Û ( −1) denotes the CS reconstruction of the spatiotemporal fMRI data at ( − 1)th FOC iteration, and ‖ ⋅ ‖ denotes the Frobenius norm. The performance of k-t FOCUSS algorithms with the proposed stopping criterion for FOC parameters was evaluated via subject-based leave-one-out cross-validation (i.e., the FOC value for data from each subject was determined based on a training set consisting of data from the other subjects). The effects of CG and were investigated separately for verification with the determined optimal 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 = ‖y −ŷ‖ 2 2 was used as a measure to observe data fitting and convergence in each k-t FOCUSS algorithm with increase in CG , where y denotes the -space time-series data of the sampled -space lines andŷ denotes the CS reconstruction of the -space time-series data for the corresponding -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 ,û denotes the CS reconstructed spatiotemporal fMRI data at time frame , denotes the total number of time frames, and denotes the total number of image pixels at each time frame.

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 0 ) according to the -statistics map of the original full-sampled data were chosen as the ROI for further analysis. New ROIs were defined for each dataset.

Quantitative Analysis.
Frame-by-frame normalized MSE, -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 was performed using the following equation: Student's -test was performed for each dataset to statistically analyze fMRI data and generate the -statistics functional map. The -score is calculated on a pixel by pixel basis over time as follows: -score = − where ⋅ denote the mean, ⋅ denote the standard deviation, and ⋅ denote the length of the baseline time-series and activation time-series , respectively. The -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, thestatistics 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.

Determination of
. FOC values determined via subject-based leave-one-out cross-validation for Algorithm 1 were 4, 3, 3, and 3 for Pattern R , Pattern G , Pattern GR , and Pattern GRC1 , respectively, and those for Algorithm 2 were 5, 4, 4, and 4 for Pattern R , Pattern G , Pattern GR , and Pattern GRC1 , respectively. Identical 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 FOC values for each subject data, to evaluate the performance of -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.

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 -space lowfrequency 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 -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., Pattern GR and Pattern GRC1 ) 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 -space low-frequency information showed increase of artifacts (Figure 2(c)).
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 Pattern R did not show any meaningful functional activations ( Figure 6(a)). The fMRI maps from downsampling with only -space low-frequency information showed significant blurring in the activation region (Figure 2(d)). The fMRI maps from Pattern GR (Figure 6(c)) and Pattern GRC1 (Figure 6(d)) were closer to the original fMRI maps than those from Pattern G (Figure 6(b)) in terms of preserving details in activation foci shift, presumably due to the inclusion of appropriate high frequency -space information.
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 Pattern R 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). Pattern GRC1 displayed the highest ROC performance than other sampling patterns including downsampling with only -space low-frequency information.
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 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).

Reconstruction with Algorithm 2: k-t FOCUSS with KLT.
Results of Algorithm 2 in terms of baseline images, frame-byframe 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).

Investigation of and Effect.
Representative results from CS reconstruction of Pattern GRC1 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 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).

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 -space center region such as Gaussianweighting or inclusion of -space center 1 line was better for CS, due to the fact that most energy is located in the -space center region and also due to incoherent aliasing effects from variable-density sampling [13,16]. The reconstruction results from Pattern GR and Pattern GRC1 also verify that a variation in sampling scheme with more -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 morespace low-frequency information implies less -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 -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 eventrelated 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 kt 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][33][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 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 -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 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 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 CG and , the only issue for the application of k-t FOCUSS algorithms for high field bSSFP fMRI data lies in the choice of FOC . The 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 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 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 Frame-by-frame normalized MSE(t) = ‖u t −û t ‖ 2 2 /‖u t ‖ 2 2 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. al. [35] discovered a simple method to suppress the eddy current effect by pairing two consecutive -space lines. By incorporating this idea, a simulation study was performed to see the effect of pairwise downsampling scheme. A paired Pattern GRC1 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.
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, respiratoryinduced noise, cardiac pulsation, and subject motion [37][38][39] and can also improve the activation detection sensitivity instatistics [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   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) Pattern R , (b) Pattern G , (c) Pattern GR , and (d) Pattern GRC1 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.
CS [18,[40][41][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 2fold 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 Frame-by-frame normalized MSE(t) = ‖u t −û t ‖ 2 2 /‖u t ‖ 2 2 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.
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.

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 -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.

Review of -FOCUSS
k-t FOCUSS with Temporal Fourier Transform. k-t FOCUSS is a recent CS algorithm developed for the reconstruction of dynamic image data [16,18]. As the name indicates, it is based on FOCal Underdetermined System Solver (FOCUSS) algorithm and utilizes random sampling in the -domain [16,43,44]. Here, the basic structure of the algorithm will be speculated. For simplicity, only the case of Cartesian -space trajectory will be discussed, with downsampling only in the phase-encoding direction and full sampling in the frequencyencoding direction.
In a discrete setup, the measurement-signal relationship at time is where F ∈ C × denotes a downsampled Fourier transform (FT) along the phase-encoding direction, y ∈ C denotes the space measurement vectors, u ∈ C denotes the image vector at time , and denotes the number of time frames. If image content varies periodically over time, we can sparsify its temporal variation using × FT matrix Φ such that the corresponding coefficients {x } become sparse: where y = vec(Y) ∈ C and x = vec(X) ∈ C . In k-t FOCUSS, the unknown signal x is further decomposed as where x 0 denotes a predicted signal (such as temporal mean) and Δx denotes residual signal that needs to be reconstructed using CS. Accordingly, the CS formulation is given by As an optimization method for (A.6), k-t FOCUSS employs weighted-2 minimization or FOCUSS algorithm by converting Δx = Wq, which provides the following unconstrained form of cost function:   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) Pattern R , (b) Pattern G , (c) Pattern GR , and (d) Pattern GRC1 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.   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 Pattern GRC1 and optimal 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 .  Here, q = [ 1 , . . . , ] denotes the estimated signal x from the previous iteration with = . Then, the optimal solution of the problem is found by calculating the following: x = x 0 + ΘA (AΘA + I) −1 (y − Ax 0 ) , (A.9) where Θ = WW . Since the direct calculation of (A.9) is computationally demanding due to the matrix inversion of a large size matrix, conjugate gradient (CG) method is used to minimize the cost function (A.7).
A generic form of the implementation of k-t FOCUSS that utilizes temporal FT Φ is summarized in Algorithm 1.
k-t FOCUSS with KLT. Even though k-t FOCUSS has been developed as above utilizing temporal FT as Φ in (A.3), other temporal transformations could also be used to efficiently sparsify the signal [16]. One example is the utilization of Karhunen-Loéve transform (KLT), which is also known as principal component analysis (PCA) [45].
KLT or PCA is a data-dependent mathematical procedure that uses an orthogonal transformation to convert possibly correlated elements of the data into a set of linearly uncorrelated components called "principal components. " The transformation leads to a result in which the first principal component accounts for the largest possible variance of the data, and each succeeding component has the next largest variance possible under the restriction that it is orthogonal (i.e., uncorrelated) to the preceding components. The transform is known to compact most of the energy into a small number of expansion coefficients [45] and thus is ideal in CS perspective [16].
The application of KLT within k-t FOCUSS requires calculation of a covariance matrix from a trained dataset. In this paper, the result from k-t FOCUSS that utilizes temporal FT using Algorithm 1 is used for the calculation of covariance matrix. Once the training dataset is defined, the covariance matrix is constructed as follows: (A.10) Then, the eigenvectors of C can be used for the KL transform; that is, After updating Φ, we can perform additional k-t FOCUSS iterations. The algorithm is summarized in Algorithm 2.