Abstract

In transition-band steady-state free precession (SSFP) functional MRI (fMRI), functional contrast originates from a bulk frequency shift induced by a deoxygenated hemoglobin concentration change in the activated brain regions. This frequency shift causes a magnitude and/or phase-signal change depending on the off-resonance distribution of a voxel in the balanced-SSFP (bSSFP) profile. However, in early low-resolution studies, only the magnitude signal activations were shown. In this paper the task-correlated phase-signal change is presented in a highresolution (1 ؋ 1 ؋ 1 mm 3 ) study. To include this phase activation in a functional analysis, a new complex domain data analysis method is proposed. The results show statistically significant phase-signal changes in a large number of voxels comparable to that of the magnitude-activated voxels. The complex-data analysis method successfully includes these phase activations in the activation map and thus provides wider coverage compared to magnitude-data analysis results. Since its first inception in the early 1990s, blood oxygen level-dependent (BOLD) functional MRI (fMRI) has had a tremendous impact on human brain research. Most studies performed a functional data analysis using the magnitude data of the time-series images while ignoring phase-signal changes. Recently a few high-resolution studies explored the phase information and revealed that the task-correlated phase-signal change can provide additional information regarding brain activation (1,2). Investigators found that a large draining vein induced phase-signal changes in addition to magnitude-signal changes by creating a relatively uniform frequency shift within a voxel. This information was utilized to determine the oxygen saturation change in veins (1) and to suppress BOLD signals from large vessels (2). However, in transition-band SSFP, or blood oxygenation sensitive steady state (BOSS) fMRI, a recently developed functional imaging method that provides greater signal-to-noise ratio (SNR) efficiency with less imaging distortion (3-6), the origin of the phase-signal change differs from that in conventional gradient-recalled echo (GRE) fMRI (5). Moreover, this task-correlated phasesignal change was always expected, since a sharp phase transition near the on-resonance frequency (also at every multiple of 1/TR) is considered to be a primary source of the functional contrast in transition-band SSFP fMRI. However, early low-resolution studies measured functional activations from the magnitude signals (5). This has raised questions regarding the existence of the phase activation and the contribution of the SSFP phase profile to the functional contrast. In this article we present the results of a 3D high-resolution (1 mm 3 ) functional study to show the presence and characteristics of task-correlated phase-signal changes in transition-band SSFP fMRI. Simulation results are also included to further elucidate the sensitivity and characteristics of the phase-signal changes. If functional contrast exists in both magnitude and phase signals, and if both signals provide spatially localized information, an analysis based only on magnitude cannot fully exploit the data. Some voxels that contain phase activation without magnitude activation will not be detected, and some voxels that contain both activations will show less activation in the magnitude-only analysis. Therefore, to acquire more reliable results, a complex domain data analysis method that encompasses both magnitude-and phase-signal activations is necessary. Recently several complex-data analysis methods (7-15), most of which are based on a generalized likelihood ratio test (GLRT) (13), have been proposed. Here we propose a new method that combines a Hotelling's T 2 -test (which can be derived from GLRT) (16) with a generalized linear model (GLM) (17) to calculate the statistical significance of the activation from complex data. This method is computationally efficient compared to the previously proposed methods, and generates full activation maps that include both magnitude and phase activations. THEORY Functional Contrast of Transition-Band SSFP fMRI Unlike the contrast mechanism of T 2 * dephasing in conventional GRE fMRI, transition-band SSFP fMRI is based on a bulk frequency shift induced by a fractional oxygen saturation change of hemoglobin where ␥ is the gyromagnetic ratio (2.678 ϫ 10 8 rad T -1 s -1 ), Hct is a fractional hematocrit in blood (0.4), ⌬ is the susceptibility difference between fully oxygenated and fully deoxygenated red blood cells (0.27 ppm from Ref. 905 20), Y is the fractional oxygen saturation of hemoglobin, R is the radius of the vessel, and r is a vector from the center axis of the vessel on an orthogonal plane of the vessel. The magnitude and phase signals of a voxel are determined by the intra-and extravascular frequency shifts, the fractional volume of the vessels, and the off-resonance frequency distribution. If a voxel contains n veins whose volumes are small (and the veins are distant from each other), and the off-resonance frequency distribution of the voxel is assumed to be uniform, then the magnetization of the voxel (M) becomes where ␣ is the total volume of all veins (including extravascular spaces), f off is the off-resonance frequency of the voxel, ␣ i is the volume of the i th vein (including the extravascular space), ⌬f i is the fractional frequency shift from the off-resonance frequency induced by the i th vein in the d␣ i space, and P is the magnetization profile of balanced SSFP (bSSFP; If a vessel is located near the resonance frequency, the sharp phase transition will create a large phase-signal change in accordance with the bulk frequency shift. As a result, a task-correlated phase-signal change is expected. Based on the orientation of the vessel, the intravascular frequency shift can cause either a positive (|| Ͻ 0.96 or 2.19 Ͻ || Ͻ 3.14 radian) or a negative (0.96 Ͻ || Ͻ 2.19 radian) phase change, resulting in two different phase correlations. The extravascular frequency shift can also contribute to the net phase accrual due to the nonlinear phase profile. However, if a voxel is large and only a fractional volume of the vessels shows oxygen concentration changes, and the vessels are randomly oriented, then the task-correlated phase-signal change will be decreased due to partial-volume effects and incoherent phase accruals. As a result, in low-resolution studies the phase-signal change is difficult to observe, even in the phase transition band. Outside of the phase transition band, the bulk frequency shift (⌬f) modulates the magnitude signal from the magnitude profile with little or no phase-signal change. Both a positive and a negative magnitude change can be detected based on the off-resonance frequency and vessel orientation. For instance, if the off-resonance frequency shift is positive and a vessel is parallel to the field, the voxel will show a positive correlation with the given stimulus, whereas a vein tilted by a /2 radian will show a negative correlation. Hotelling's T 2 -Test In conventional fMRI analysis, where there is only one variate (magnitude data), a Student's t-test can be used to determine whether the null hypothesis is rejected (i.e., the voxel shows activation) or is not rejected (i.e., the voxel shows no activation) for a given threshold. However, in the complex data analysis, the number of variates becomes two (i.e., real and imaginary or magnitude and phase); hence, one must generalize the Student's t-test to a bivariate test to evaluate the null hypothesis in the complex domain. This generalization can be achieved using a Hotelling's T 2 -test (16) that compares the sample mean vector with the expected mean vector, and performs the null hypothesis test on a multivariate domain as described below. Similarly to the Student's t-test statistic t ϭ ͑x Ϫ ͒/͑ / ͱn͒, (where n is the number of samples, x is the sample mean, is the hypothetical population mean, and is the sample standard deviation (SD)), the Hotelling's T 2 -test is defined as where n is the number of samples, x is the sample mean vector, is the hypothetical population mean vector, and Ŝ is the sample covariance matrix (an unbiased estimator). When the sample has two variates (for example, real and imaginary), x and become 2 ϫ 1 vectors, whereas Ŝ becomes a 2 ϫ 2 matrix. If all of the variates (in this case, both the real and imaginary parts of the data) are distributed by the Gaussian distribution, the T 2 statistic will follow an F-distribution with m and n-m degrees of freedom (where m is the number of variates in the sample (i.e., m ϭ 2 in the complex data) and n is the number of samples). This T 2 statistic can [5] A P-value can be obtained from the cumulative distribution function (CDF) of F m,n-m evaluated at T 2 ͑n Ϫ m͒/͓͑n Ϫ 1͒m͔. Complex Data Analysis Method To incorporate the hemodynamic response function into this T 2 -test, a GLM (17) is utilized to estimate the mean vectors of the activation and baseline states, as well as the covariance matrix. The complex time-series data are decomposed into real and imaginary axes. To find the mean vector of each state, one structures the design matrix (X) of the GLM by a constant vector (1 ϭ [1 1 . . .. 1] T , a real n ϫ 1 vector) and a reference waveform vector (h, a real n ϫ 1 vector, the convolution of a stimulus pattern and a hemodynamic response function). Equation where y ϭ y r ϩ iy i (a complex n ϫ 1 vector) is the time-series data of one voxel, ␤ r and ␤ i (a real 2 ϫ 1 vector each) are the parameters of GLM, and r and i are residual errors (a real n ϫ 1 vector each). One can easily incorporate other terms, such as a linear drift, by adding more vectors and parameters into the model (see Appendix B). The least-square estimates of ␤ r and ␤ i are (X T X) Ϫ1 X T y r and (X T X) Ϫ1 X T y i , respectively. The activation level or the mean difference between the two states (x Ϫ ) becomes ͓␤ r2 ␤ i2 ͔ T , which are the parameters of the reference waveform vector in real and imaginary axes. Setting the contrast (v) to the reference waveform vector (i.e., , we calculate the covariance matrix of the contrast from the residual errors ( r and i ) and the design matrix as follows: [8] This closely resembles the variance of a contrast in a 1D GLM test whose value is 2 v With these values calculated from the GLM method, the T 2 value can be obtained as follows: The P-value is found from Eq. [6] with m ϭ 2. MATERIALS AND METHODS Simulated Data Analysis Using the Complex Data Analysis Method Four different data sets of the complex data (no contrast, contrast in magnitude, contrast in both magnitude and phase, and contrast in phase) were generated to validate the proposed complex data analysis method. Each complex data set was generated by Eq. [10] with the proper SNR and contrast-to-noise ratio (CNR) values: where y (a complex 50 ϫ 1 vector) is a voxel time-series with 50 points, X is the same design matrix as in Eq. [7], the reference waveform vector (h) of X is a boxcar block design with a duration of 10 samples, n r and n i are real 50 ϫ 1 vectors distributed by N(0,1), and ␤ r ϭ ͓␤ r1 ␤ r2 ͔ and ␤ i ϭ ͓␤ i1 ␤ i2 ͔ are the parameters in the real and imaginary axes, respectively. Since the noise is normalized, the parameters ␤ r1 and ␤ i1 represent the SNRs in the real and imaginary axes, whereas ␤ r2 and ␤ i2 represent the CNRs in the real and imaginary axes. The absolute values of the complex data were used to generate the magnitude data. We analyzed the magnitude data using the GLM method with the same design matrix (X) to compare the results. A total of 100,000 voxels were simulated in each data set. The SNRs were set to 10 (i.e., ␤ r1 ϭ ␤ i1 ϭ 10). The first data set was Gaussian noise data with ␤ r2 ϭ ␤ i2 ϭ 0. This test was intended to show that the null hypothesis is not rejected when no activation exists. The second data set (␤ r2 ϭ ␤ i2 , both simulated at increasing values from 0 to 1.4) was generated to simulate the case in which the contrast exists only in magnitude. In the third data set the same contrast existed both in magnitude and phase (␤ r2 ϭ 0 whereas ␤ i2 increased from 0 to 1.4 ͱ2). For the last data set, ␤ r2 was changed from 0 to 1.4, while ␤ i2 was kept as the opposite sign of ␤ r2 to simulate the activation approximately in phase with the same amount of the complex contrast as in the previous two data sets. Experimental Data Acquisition All experiments were performed on a 1.5 T GE EXCITE system (40 mT/m and 150 mT/m/ms) with a 3-inch receive-only surface coil, except for a high-resolution reference scan that was acquired using an eight-channel head coil (MRI Devices Corp., USA). Five subjects, who provided written consent (approved by Stanford University), were immobilized by pads and instructed to avoid any voluntary motions. A high-resolution reference scan (spoiled gradient-echo (SPGR), FOV ϭ 22 cm, resolution ϭ 1 ϫ 1 ϫ 1 mm 3 , TR ϭ 11.7 ms, TE ϭ 5.1 ms, flip angle ϭ 25°, number of excitations (NEX) ϭ 3) that covered the entire brain was acquired using the head coil on a separate day. In the main experiments using the surface coil, an axial slab (2-cm thickness) of the lower occipital lobe was targeted after the brain was localized with a 3D localizer sequence. An intermediate-stage anatomical scan (SPGR, FOV ϭ 16 cm, resolution ϭ 1 ϫ 1 ϫ 1 mm 3 , TR ϭ 12 ms, TE ϭ 4.2 ms, flip angle ϭ 25°, NEX ϭ 6) that helped to realign the high-resolution reference scan to the functional Complex Analysis in SSFP fMRI 907 data was performed at the same resolution and FOV as in the functional scans. To identify the locations of large vessels, a 3D spiral trajectory high-resolution venogram (GRE, FOV ϭ 16 cm, resolution ϭ 0.5 ϫ 0.5 ϫ 1 mm 3 , TR ϭ 70 ms, TE ϭ 40 ms, flip angle ϭ 25°, NEX ϭ 12) was obtained by using a flow-compensated venogram (21) sequence. Linear shimming was targeted at the occipital lobe of the brain by a custom targeted shim program before the functional scans. For the functional studies, a 3D stack-ofspirals sequence (bSSFP, FOV ϭ 16 cm, resolution ϭ 1 ϫ 1 ϫ 1 mm 3 , TR ϭ 15 ms, TE ϭ 1.5 ms, flip angle ϭ 5°, number of interleaves ϭ 10, 18 slices for subjects 1-3 and 16 slices for the others) was utilized to cover a 3D volume every 3 s Experimental Data Analysis All functional data were inspected for subject motion. The third scan (⌬f center ϭ 3 Hz) of subject 4 and the second and third scans (⌬f center ϭ -3 Hz and 3 Hz) of subject 5 showed significant movements (Ͼ1 mm) and were therefore discarded from further analysis. A small shift (one voxel) between the functional scans was found in subject 3, which was fixed by aligning the data to the first functional scan. The intermediate-stage anatomical images were also aligned to the first functional scan. The signal attenuation in this intermediate-stage anatomical image, induced by the surface coil sensitivity, was removed with the use of a custom program. The high-resolution reference scan was then aligned to the intermediate-stage anatomical data using SPM5 (22). The alignment was further refined manually. The gray matter (GM) regions were identified from the reference scan. The venogram result was aligned to the high-resolution reference scan. We identified the locations of large veins from both the magnitude and phase venogram results by carefully tracking the vessel geometries from all of the slices. To display the results, the reference scan and z-score maps were enlarged to match the resolution with the venogram. The magnitude and phase data were generated from the complex functional data. The phase data were unwrapped in the time-series for each voxel and 2 was added to make positive-signed data. After that, both data sets were processed individually to create "magnitude-only" and "phase-only" z-statistics maps using FEAT FSL (23). For these null-hypothesis tests, the noise should be assumed to be Gaussian. This assumption becomes valid in both magnitude and phase when the SNR is high (8). With the exception of high-pass filtering (cutoff frequency ϭ 0.022 Hz), no other prestatistical process was performed. Unthresholded z-statistics maps were generated from FEAT based on each voxel. Since the negative correlations were meaningful, the thresholds were set on both the positive and negative sides of the z-score distribution with a two-tailed P-value of 0.01. After thresholding, isolated activations in the positive and negative z-statistics maps were removed separately. Finally, the absolute values of the z-statistics maps were used to generate the activation maps. These activation maps were further masked by the mask generated in the magnitude-only data (10% thresholding, FEAT) with the skull area removed. The singlefrequency acquisition analysis was performed on the data whose center frequency shift (⌬f center ) was equal to 0 Hz. In the multifrequency acquisition analysis that combined the three different center frequency results using the maximum z-score projection method (5), the threshold was increased to P Ͻ 0.005. For the complex data analysis, the functional data were decomposed into real and imaginary time-series. Slow signal drift was then removed by the same filter that was used previously. In each voxel the T 2 -value was calculated from Eq. [9] as described in the Theory section. The reference waveform vector for the design matrix was the convolution of the stimulus pattern and a delayed-Gamma function that was generated by FEAT (the same waveform vector as that in the magnitude and phase analyses). The P-value was calculated in each voxel from the F-distribution (Eq. [6]), and the activation maps were generated by thresholding with a one-tailed P-value of 0.01 for the single-frequency acquisition analysis. After thresholding, we converted the P-values into z-scores by finding the equal probability in the one-tailed Gaussian distribution N(0,1) for the color-coded activation maps. Isolated activations were removed from the activation maps and the results were masked using the same mask used in the magnitude and phase activation maps. After calculating the activation maps for each run, we combined the three different center frequency scans using the maximum z-score projection. The same threshold (P Ͻ 0.005) that was used in the magnitude and phase multifrequency acquisition analysis was utilized to threshold the activation maps before the projection. RESULTS The simulation results are shown in The experimental results reported in the following paragraphs refer to subject 1 (the results from all five subjects are summarized in The maximum z-score projected phase activation maps (P Ͻ 0.005) from the three scans reveal high z-score voxels near the veins (highlighted in yellow) that can be identified from the venogram (blue circles in In The maximum z-score projected activation maps from the complex-data analysis are shown in DISCUSSION Task-Correlated Phase-Signal Change In this study we observed the expected task-correlated phase-signal change in a transition-band SSFP fMRI experiment using a 3D high-resolution (1 ϫ 1 ϫ 1 mm 3 ) acquisition. A significant difference compared to previous lowresolution studies is that large task-correlated signal changes were observed in both magnitude and phase. The magnitude results were similar in range to those reported in Ref. 6; however, the large phase-signal change and the large number of voxels that showed significant phasesignal changes were observed for the first time, to our knowledge, in the transition-band SSFP fMRI experiment. To further validate these signal change levels, we performed a simulation using Eqs. [1]- More simulation results to elucidate the nature of the phase-signal change are shown in Complex Analysis in SSFP fMRI 911 ent angles (0 Յ Յ , where is an angle of the vein with respect to a B 0 field). All of the sequence parameters were the same as before (TE ϭ TR/2). The intravascular phase signal decreased as the angle increased from 0 to 0.96 (rad) and became negative from 0.96 to /2 (rad). In the extravascular case, the phase signal increased (both positively and negatively depending on the off-resonance frequency) as the angle increased from 0 to /2. From /2 to , both intra-and extravascular phase activations reversed the patterns from 0 to /2 (see [1]), the saturation of the phase activation will start at a higher field strength in other vessel orientations. Moreover, the contributions from the extravascular frequency shift will also have different saturation patterns. In the capillary case, the oxygen concentration change (Y) is smaller than that of the vein, resulting in smaller intravascular frequency shifts that in turn reduce the saturation effect. One criterion for including the phase-signal change in the activation map is that it must show functionally localized signals. In a study by Menon (2), the source of the phase-signal change was assumed to originate primarily from the large veins. As a result, the phase-signal change was not included, but rather was utilized to suppress the large vein signals. In his experiment the data were collected at 4T, where (theoretically) the bulk frequency shift is 2.7 times greater than it is at 1.5T. However, the phase changes in Ref. 2 (0.085 radian in visible veins and 0.028 radian in other activated voxels) appear to be smaller compared to our results (a maximum 0.482 radian phase change with 0.107 radian, on average, in the single-frequency analysis), presumably because the phase change is much greater in transition-band SSFP fMRI due to the sharp phase transition band, and also because we employed a smaller voxel size. As also mentioned in Ref. 2, the minimum detectable phase change depends on the SNR. In our experiment the SNR measured in the nonactivated voxels (|z| Ͻ 0.67) was 12.3 within the mask (roughly covering the brainstem to the back of the brain) with a single-pixel, temporalphase SD of 0.114. Hence, the minimum detectable phase change in our experimental setting (P Ͻ 0.01 and n ϭ 45) was 0.041 radian. This minimum threshold agrees with the histogram result of the signal changes from the phase activation voxels whose minimum detected signal level was 0.031 radian 912 Lee et al. phase-activated voxels was comparable to the number of magnitude-activated voxels (which can be interpreted to mean that the sensitivity of the phase activation is similar to the sensitivity of the magnitude), indicate that for the proposed method it may be beneficial to include the phase-signal change in the activation maps. There are several ways to improve the minimum detection threshold. First, one can increase the total acquisition time to increase the statistical power. Since our experiment was relatively short (2 min 15 s), increasing the total scan time by a factor of 4 (9 min) is still a reasonable approach. Another method is to use a higher-field-strength system. The bulk frequency shift of vessels increases linearly with the field strength, resulting in increased phasesignal changes. In addition to this increased contrast, the SNR also increases linearly, resulting in a quadratic increase of the detectability. A third method is using a realtime, respiration-compensation technique. Lee et al. (27) demonstrated a significant increase in the z-score by reducing the respiration-induced the SSFP profile shift. This respiration compensation is especially effective for increasing the contrast and reducing the noise interference because the respiration-induced B 0 field modulates the SSFP transition profile, causing a large signal interference and time-varying functional contrasts in transition-band SSFP fMRI. By combining these methods, one can obtain much more localized activation information from the phase signal, which makes the proposed complex data analysis method more beneficial for transition-band SSFP fMRI. Since both phase and magnitude signal changes are greater in large veins, both the complex-and magnitudedata analyses will show greater signal changes in the veins. However, the signal changes from large veins are usually of no interest because they impede the localization ability of high-resolution studies. Therefore, it is desirable to remove these large signal changes from the activation maps. This can be done with the use of diffusion-weighted suppression techniques that remove intravascular signals using bipolar diffusion gradients Once the large vein signals are removed, the complex-data analysis method will provide greater benefits by including more localized brain activations in the analysis. Complex Data Analysis Method When written in polar coordinates, the magnitude and phase contrasts are orthogonal to each other; therefore, a univariate test designed for testing changes in magnitude, by definition, will be insensitive to changes in phase. A change that occurs in magnitude and phase simultaneously can only be captured fully by a complex bivariate test. On the other hand, if the effect truly exists only in the magnitude, the magnitude-only test will have more power. This is because in order to prevent false positives in the direction of the phase, the complex test allocates a null probability mass in that direction, which is wasted if the change occurs only in the magnitude. As a result, the magnitude-only test yields a relatively lower threshold (and therefore higher power) for the same significance level. This explains the simulation described in the Results section, in which the complex-data analysis outperformed the magnitude data analysis except when the contrast existed strictly in magnitude. A mathematical analysis of the performance of the two methods is presented in Appendix A. The performance of the analysis methods was measured by the required SNR to detect a certain contrast level at a given probability of detection (P D ) and probability of false alarm (P FA ). [A1] and [A2]. The dark red straight band in [A3]). The required SNR for the magnitude data analysis is higher than that of the complex data analysis for most of the contrasts except for the dark blue areas (29% of the total plane) between the dashed white lines, where the contrasts are primarily aligned in the magnitude direction. This result indicates the superiority of the complex data analysis for most of the contrasts, which was also observed in the simulation and experimental results The proposed complex data analysis can be intuitively understood as follows: If the noise is assumed to be distributed as a 2D Gaussian with N(0, 2 I), where I is a 2 ϫ 2 identity matrix, Ŝ Ϫ1 (in Eq. [4]) becomes I/ 2 . In this case the T 2 value becomes proportional to ͑x Ϫ ͒ T ͑x Ϫ ͒, which can be seen as the distance between the sample mean vector and the population mean vector in a complex plane. In other words, if one uses a simple boxcar model, ignoring the transients of the hemodynamic response function, the signal level difference (or contrast) in fMRI data can be defined as the distance between the mean vector of the activation state and the mean vector of the baseline in the complex plane. In Appendix B a theoretical proof is given to demonstrate that our complex data analysis method is equivalent to the GLRT method (9) when the design matrices are the same for the magnitude and phase. One of the advantages of our method is that it is computationally efficient. Since it is based on the T 2 -test with the GLM method, the computational time for a sample data took only 2.4 times longer compared to the magnitude data analysis. On the other hand, the GLRT method requires iterations in each voxel to search for the optimum parameter values, which can potentially lead to a non-negligible increase in computational time for fMRI data. Future Work In addition to the bulk frequency shift, the BOLD effect includes a T 2 change that originates from the protons diffused by red blood cells (intravascular) and the microgradients between the vessel and the surrounding tissue (extravascular). In SSFP-based fMRI, this T 2 change can modulate the magnetization profile of bSSFP. A voxel with large veins can experience a T 2 on the order of a 10% change primarily from the intravascular effect, whereas the T 2 change in the GM is negligible due to the small volume of blood at 1. 5T (25,32). At higher field strengths, however, the increased extravascular effect can make the T 2 change significant in the GM. Therefore, it is important to consider the profile modulation from the T 2 change as a source of the functional contrast at higher field strengths. Another possible application of the proposed complex data analysis is conventional GRE-based fMRI in high-field experiments. At high field strengths (Ն 7T), where T 2 * shortening is significant in large veins (29), the phasesignal change from the GM will show a higher signal correlation compared to that obtained at a lower field strength. In this case the proposed complex-data analysis method will be beneficial for detecting signal changes from both the magnitude and phase. CONCLUSIONS In this paper we have presented the results of an isotropic 1-mm resolution 3D (160 ϫ 160 ϫ 18 mm 3 ) transition-band SSFP fMRI experiment. These high-resolution results provide both magnitude and phase data, depending on the off-resonance frequency and vessel properties of a voxel. To include both functional contrasts, a new complex data analysis method based on the T 2 -test is proposed. This method includes both magnitude and phase activations to compensate for the missing phase activations in the magnitude data analysis method. The method will provide more benefits at higher field strengths after large-vessel signals are removed. It can also be adopted in conventional fMRI analysis to include the task-correlated phase-signal change in high-resolution or high-field experiments. ACKNOWLEDGMENTS