Abstract
Conventional functional magnetic resonance imaging (fMRI) technique known as gradientrecalled echo (GRE) echoplanar imaging (EPI) is sensitive to image distortion and degradation caused by local magnetic field inhomogeneity at high magnetic fields. NonEPI sequences such as spoiled gradient echo and balanced steadystate free precession (bSSFP) have been proposed as an alternative highresolution fMRI technique; however, the temporal resolution of these sequences is lower than the typically used GREEPI fMRI. One potential approach to improve the temporal resolution is to use compressed sensing (CS). In this study, we tested the feasibility of kt FOCUSS—one of the high performance CS algorithms for dynamic MRI—for nonEPI fMRI at 9.4T using the model of rat somatosensory stimulation. To optimize the performance of CS reconstruction, different sampling patterns and kt FOCUSS variations were investigated. Experimental results show that an optimized kt FOCUSS algorithm with acceleration by a factor of 4 works well for nonEPI fMRI at high field under various statistical criteria, which confirms that a combination of CS and a nonEPI sequence may be a good solution for highresolution 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 leveldependent (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 gradientrecalled echo (GRE) echoplanar 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. NonEPI sequences such as spoiled gradient echo or balanced steadystate 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 GREEPI.
One solution to overcome the low temporal resolution of nonEPI sequences is to adopt parallel imaging technique [10–12]. Although proven useful, the usage of parallel imaging results in reduced signaltonoise 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 nonEPI sequences is to use compressed sensing (CS) [13–15]. CS theory states that it is possible to reconstruct an aliasingfree 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 kt FOCUSS by employing random sampling pattern in t space and by using various sparsifying temporal transforms such as Fourier transform (FT) and KarhunenLoé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 GREEPI fMRI: ordinary GREEPI [19, 20] and spiral scan GREEPI [21]. Despite its application, GREEPI 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 GREfMRI may increase statistical performance of activation detection [22]. NonEPI sequences such as bSSFP or GRE may work better with CS algorithms than GREEPI, since the sequences utilize different RF excitations for each TR.
The signal dynamics of steadystate 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 nonEPI fMRI at 9.4T using the model of rat somatosensory stimulation. Fully sampled data with highresolution spoiled gradient echo and 4 independent passband bSSFP fMRI each with different phasecycling angles (, , , and ) underwent retrospective downsampling and reconstruction using kt 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 fullsampled data. The potential for improving the temporal resolution of nonEPI 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 SpragueDawley 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 (RSP1002, 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 passband bSSFP studies were performed with TR/TE = 20/10 ms and 10/5 ms, respectively. The bSSFP fMRI studies were performed with four different phasecycling 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.
2.2. 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., spacetemporal) 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 fullsampled space dataset for retrospective downsampling before the CS reconstruction procedure.
(a)
(b)
(c)
(d)
(e)
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 lines were sampled along the phaseencoding (PE) direction to maintain the downsampling factor of 4 (where 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 , with (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 number of PE lines according to the above Gaussian probability distribution and subsequently sampling number of PE lines according to a uniform random probability distribution ( 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., ) 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 TimeSeries MRI Data.
For downsampling factor (), one has the following.(A)Uniform Random Sampling Mask(1)Generate uniform probability distribution along the PE dimension of space data.(2)Sample number of PE lines according to the uniform probability distribution.(3)Repeat steps (1)(2) for each time frame.(B)Gaussian Random Sampling Mask(1)Define in relation to .(2)Generate Gaussian probability distribution along the PE dimension of space data.(3)Sample number of PE lines according to the Gaussian probability distribution.(4)Repeat steps (1)–(3) for each time frame.(C)Mixture of Gaussian and Uniform Random Sampling Mask (Gaussian : Random = : )(1)Define in relation to .(2)Generate Gaussian probability distribution along the PE dimension of space data.(3)Sample number of PE lines according to the Gaussian probability distribution.(4)Generate uniform probability distribution along the PE dimension for the remaining PE lines.(5)Sample number of PE lines according to the uniform probability distribution.(6)Repeat steps (1)–(5) for each time frame.(D)Mixture of Gaussian and Uniform Random Sampling Mask with Full Sampling of Center 1 Line (Gaussian : Random = : )(1)Sample PE line from the space center (e.g., ).(2)Repeat all steps from (C) with number of PE lines.
The Gaussian probability distribution was considered as a derivative form of space centerweighted downsampling pattern, with stronger weighting on space lowfrequency 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 fullsampled 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 , , , and , respectively.
2.3. CS Algorithms
Two variations of kt 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 kt 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 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 was used to update once.


2.4. kt FOCUSS Parameters
For both algorithms, weighting matrix power factor () of was used to find the sparse solution equivalent to the solution of CS [16]. Conjugate gradient (CG) iteration number () of was considered sufficient and was used for both Algorithm 1 () and Algorithm 2 (), based on previous application with kt FOCUSS [25]. Regularization factor () of was used for Algorithm 1 () [16] and was used for Algorithm 2 () [22].
Both Algorithms 1 and 2 were optimized based on variation in the FOCUSS iteration number () parameters (i.e., and , resp.). The following stopping criterion is used to determine the optimal value from each dataset in training phase:where denotes the CS reconstruction of the spatiotemporal fMRI data at FOC iteration, denotes the CS reconstruction of the spatiotemporal fMRI data at FOC iteration, and denotes the Frobenius norm. The performance of kt FOCUSS algorithms with the proposed stopping criterion for parameters was evaluated via subjectbased leaveoneout crossvalidation (i.e., the value for data from each subject was determined based on a training set consisting of data from the other subjects).
The effects of and were investigated separately for verification with the determined optimal value (i.e., value found during the training phase) for each subject data. The reconstruction of phasecycled 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 kt FOCUSS algorithm with increase in , where denotes the space timeseries data of the sampled space lines and denotes the CS reconstruction of the space timeseries data for the corresponding space lines. Average mean square error (MSE) of the whole timeseries data was used as a measure to observe the noise level in the reconstructed image with variation in , which was calculated as follows:where denotes the original fullsampled 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.
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 ) according to the statistics map of the original fullsampled data were chosen as the ROI for further analysis. New ROIs were defined for each dataset.
2.6. Quantitative Analysis
Framebyframe 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° phasecycling) displaying a fairly uniform signal contrast without any significant banding artifacts and showed clear activation foci in the fullsampled data.
The framebyframe 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:where denote the mean, denote the standard deviation, and denote the length of the baseline timeseries and activation timeseries , 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 signaldetection accuracy [26]. For each dataset, the statistics map generated from the original fullsampled 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 relates to specificity.
3. Results
3.1. Determination of
values determined via subjectbased leaveoneout crossvalidation for Algorithm 1 were 4, 3, 3, and 3 for , , , and , respectively, and those for Algorithm 2 were 5, 4, 4, and 4 for , , , and , respectively. Identical 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 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 phasecycling angles in the fMRI maps of fullsampled 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 phasecycled 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.
3.3. Reconstruction with Algorithm 1: kt 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 framebyframe normalized MSE values from all the CS sampling patterns were significantly lower than those from downsampling with only space lowfrequency information (Figure 5), indicating high reconstruction performance of Algorithm 1. In particular, the mixture of Gaussian and uniform random sampling scheme (i.e., and ) showed the lowest framebyframe normalized MSE values across all time frames. Overall, all Gaussianweighted 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 phasecycled bSSFP data, while downsampling with only space lowfrequency information showed increase of artifacts (Figure 2(c)).
The fMRI maps were also reconstructed well from all Gaussianweighted sampling patterns using Algorithm 1 (Figures 6(b), 6(c), and 6(d)), while those from did not show any meaningful functional activations (Figure 6(a)). The fMRI maps from downsampling with only space lowfrequency information showed significant blurring in the activation region (Figure 2(d)). The fMRI maps from (Figure 6(c)) and (Figure 6(d)) were closer to the original fMRI maps than those from (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 Gaussianweighted sampling patterns (Figures 7(b), 7(c), and 7(d)), while those from differed from the original with significantly increased temporal fluctuation (Figure 7(a)). Mean ROI time courses of images from all Gaussianweighted sampling patterns resembled those of the original data; even with only 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 Gaussianweighted sampling patterns using Algorithm 1 (Table 1). displayed the highest ROC performance than other sampling patterns including downsampling with only space lowfrequency information.
(a)
(b)
(c)
(d)
The reconstruction times of the kt 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 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).
3.4. Reconstruction with Algorithm 2: kt FOCUSS with KLT
Results of Algorithm 2 in terms of baseline images, framebyframe 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).
(a)
(b)
(c)
(d)
3.5. Investigation of and Effect
Representative results from CS reconstruction of 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 value of (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 kt FOCUSS algorithms was relatively insensitive to variation and minimal error was achieved with small values of (e.g., less than 1).
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
(d)
(e)
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 nonEPI 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 variabledensity sampling [13, 16]. The reconstruction results from and 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 more space lowfrequency 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 nonEPI 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 nonEPI fMRI.
The two algorithms of kt 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, kt FOCUSS with temporal FT performed slightly better than kt 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 spatialtemporal frequency domain may have been exploited better for such data type. Since KLT is a datadriven transform, kt 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 nonEPI 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 crossvalidation [32–34]. There are multiple parameters involved for the case of kt FOCUSS algorithm, which requires hyperparameter optimization and thus increases the computation burden [16]. Therefore, the effect of two different kt FOCUSS parameters, and , 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 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 and estimation ) as number of CG iterations increases (note that the signalmeasurement relationship is defined in (A.4) of the Appendix and can be used to find from ). Based on the results from Figures 12 and 14, fixation of 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 kt FOCUSS algorithms on high field nonEPI fMRI studies. These results agree with previous applications of kt FOCUSS where a sufficient value of and a small value of are used and are proven to perform well in highquality 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 and , the only issue for the application of kt FOCUSS algorithms for high field bSSFP fMRI data lies in the choice of . The 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 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 tradeoff between minimization of error and increase in postprocessing time must be considered appropriately before choosing the value for future studies.
Eddy currents can cause problems in bSSFP imaging with nonlinear phaseencoding 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 space lines. By incorporating this idea, a simulation study was performed to see the effect of pairwise downsampling scheme. A paired was generated with a downsampling factor of 4. The downsampling scheme and reconstructed fMRI maps from each kt 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 kt 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 tradeoff between the suppression of eddy current effect and fMRI sensitivity as well as specificity.
Application of CS to high field nonEPI fMRI can be meaningful for highresolution fMRI studies, since conventional GREEPI fMRI is sensitive to image distortion and degradation caused by local magnetic field inhomogeneity at high magnetic fields. Although the temporal resolution of nonEPI sequences is lower than the typically used GREEPI, 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 weightednorm process of the CS algorithm can reduce artifacts from scannerrelated drifts, respiratoryinduced noise, cardiac pulsation, and subject motion [37–39] and can also improve the activation detection sensitivity in 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 kt 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 nonEPI 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 tradeoff 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 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 nonEPI 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 4fold 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 nonEPI fMRI studies.
In this paper, the effect of CS is investigated through retrospective downsampling of fullsampled 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 kt 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 nonEPI fMRI data. The two kt 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 nonEPI fMRI studies. The application of CS to fMRI has great potential in practice for improvement of temporal resolution and/or spatial coverage.
Appendix
Review of FOCUSS
kt FOCUSS with Temporal Fourier Transform. kt 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 phaseencoding direction and full sampling in the frequencyencoding direction.
In a discrete setup, the measurementsignal relationship at time iswhere denotes a downsampled Fourier transform (FT) along the phaseencoding direction, denotes the space measurement vectors, 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 become sparse:where .
Then, the measurementsignal relationship becomesUsing the property , where denotes columnwise vectorization operation, we havewhere and . In kt FOCUSS, the unknown signal is further decomposed aswhere denotes a predicted signal (such as temporal mean) and denotes residual signal that needs to be reconstructed using CS. Accordingly, the CS formulation is given byAs an optimization method for (A.6), kt FOCUSS employs weighted minimization or FOCUSS algorithm by converting , which provides the following unconstrained form of cost function:whereHere, denotes the estimated signal from the previous iteration with . Then, the optimal solution of the problem is found by calculating the following:where . 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 kt FOCUSS that utilizes temporal FT is summarized in Algorithm 1.
kt FOCUSS with KLT. Even though kt 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 KarhunenLoéve transform (KLT), which is also known as principal component analysis (PCA) [45].
KLT or PCA is a datadependent 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 kt FOCUSS requires calculation of a covariance matrix from a trained dataset. In this paper, the result from kt 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: Then, the eigenvectors of can be used for the KL transform; that is, After updating , we can perform additional kt FOCUSS iterations. The algorithm is summarized in Algorithm 2.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was supported by the National Research Foundation (NRF) of Korea under project number NRF2013M3A9B2076548 and Samsung Research Funding Center of Samsung Electronics under project number SRFCIT140105.