Single-shot compensation of image distortions and BOLD contrast optimization using multi-echo EPI for real-time fMRI

www.elsevier.com/locate/ynimg NeuroImage 24 (2005) 1068 – 1079

Single-shot compensation of image distortions and BOLD contrast optimization using multi-echo EPI for real-time fMRI Nikolaus Weiskopf,a,b,c,* Uwe Klose,b Niels Birbaumer,a,d and Klaus Mathiake a

Institute of Medical Psychology and Behavioral Neurobiology, University of Tu¨bingen, Tu¨bingen, Germany Section of Experimental MR of the CNS, Department of Neuroradiology, University of Tu¨bingen, Tu¨bingen, Germany c Wellcome Department of Imaging Neuroscience, Institute of Neurology, University College London, London, UK d Center for Cognitive Neuroscience, University of Trento, Trento, Italy e Center for Neurology, Hertie Institute for Clinical Brain Research, University of Tu¨bingen, Tu¨bingen, Germany b

Received 24 April 2004; revised 17 August 2004; accepted 8 October 2004 Available online 8 December 2004

Functional magnetic resonance imaging (fMRI) is most commonly based on echo-planar imaging (EPI). With higher field strengths, gradient performance, and computational power, real-time fMRI has become feasible; that is, brain activation can be monitored during the ongoing scan. However, EPI suffers from geometric distortions due to inhomogeneities of the magnetic field, especially close to air–tissue interfaces. Thus, functional activations might be mislocalized and assigned to the wrong anatomical structures. Several techniques have been reported which reduce geometric distortions, for example, mapping of the static magnetic field B0 or the point spread function for all voxels. Yet these techniques require additional reference scans and in some cases extensive computational time. Moreover, only static field inhomogeneities can be corrected, because the correction is based on a static reference scan. We present an approach which allows for simultaneous acquisition and distortion correction of a functional image without a reference scan. The technique is based on a modified multi-echo EPI data acquisition scheme using a phase-encoding (PE) gradient with alternating polarity. The images exhibit opposite distortions due to the inverted PE gradient. After adjusting the contrast of the images acquired at different echo times, this information is used for the distortion correction. We present the theory, implementation, and applications of this single-shot distortion correction. Significant reduction in geometric distortion is shown both for phantom images and human fMRI data. Moreover, sensitivity to the blood oxygen level-dependent (BOLD) effect is increased by weighted summation of the undistorted images. D 2004 Elsevier Inc. All rights reserved. Keywords: Distortion correction; Echo-planar imaging; EPI; Functional magnetic resonance imaging, fMRI; Real-time; Single-shot; BOLD; Sensitivity; Multi-echo; Multi-image

* Corresponding author. Wellcome Department of Imaging Neuroscience, Institute of Neurology, University College London, 12 Queen Square, London WC1N 3BG, UK. Fax: +44 20 7813 1420. E-mail address: [email protected] (N. Weiskopf). Available online on ScienceDirect (www.sciencedirect.com). 1053-8119/$ - see front matter D 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2004.10.012

Introduction Echo-planar imaging (EPI) is widely used for functional magnetic resonance imaging (fMRI) of the human brain. This single-shot imaging technique offers a high temporal resolution and is relatively insensitive to motion and physiological artifacts, because a single slice is encoded in under 100 ms. It even allows for imaging of human brain function in real time (e.g., Posse et al., 2001) and neurofeedback of fMRI signals (deCharms et al., 2004; Posse et al., 2003a; Weiskopf et al., 2003, 2004; Yoo and Jolesz, 2002). In addition to fast data acquisition and processing, real-time fMRI applications such as functional localizers (Goodyear et al., 1997; Weiskopf et al., 2004) or neurofeedback demand high image quality and superior functional contrast-to-noise ratio (CNR) in order to increase sensitivity to brain activations. Single-shot multi-echo EPI (also referred to as single-shot multi-image EPI) acquires additional images as compared to conventional EPI and increases CNR (Posse et al., 1999; Speck and Hennig, 1998). Moreover, this technique can be implemented for real-time fMRI (Posse et al., 1999, 2001). Both conventional EPI and multi-echo EPI suffer from geometric distortions, if there are local differences in the resonance frequency. Distortions are most prominent along the phase-encoding (PE) direction, because it is encoded with a low effective spectral bandwidth. Distortions along other directions might occur as well but are much smaller and were disregarded in this study (Morgan et al., 2004; Munger et al., 2000). The differences in the resonance frequency, that is, frequency offsets, might result from imperfections in the magnetic field gradients used for imaging, for example, those caused by eddy currents, and gradients in the concomitant fields (Jezzard and Clare, 1999). More importantly, frequency offsets are caused by inhomogeneities of the static magnetic field which are mainly due to differences in the magnetic susceptibilities of the imaged object, for example, different susceptibility of air, bone, and brain tissue (Jezzard and Clare, 1999). Frequency offsets caused by susceptibility differences increase with the strength of the static magnetic field of the magnetic resonance (MR) scanner. Therefore,

N. Weiskopf et al. / NeuroImage 24 (2005) 1068–1079

geometric distortions are more prominent at high-field systems (z3 T), using the same readout as at low fields. These systems, however, provide an improved signal-to-noise ratio (SNR) and increased sensitivity to the blood oxygen level-dependent (BOLD) effect used in fMRI (Krasnow et al., 2003). Thus, they are used more and more frequently. Since the localization of brain activity is based on EPI data in fMRI experiments, geometric distortions lead to mislocalizations and make it difficult to assign activations to their respective anatomical location (Studholme et al., 2000). Even worse, effects such as respiration (Glover et al., 2000) and head motion (Jezzard and Clare, 1999; Ward et al., 2002) may cause dynamic changes in this mislocalization during the experiment. We aimed at the development of a real-time distortion correction which accounts for dynamic changes in the distortions. Ideally, such a distortion correction technique should be fast and robust, and it should retain or even increase BOLD CNR. Optimization of the technical CNR is crucial for fMRI, since it is the chief determinant of sensitivity to functional activations and the statistical power of the procedure. Besides advanced shimming techniques which improve global field homogeneity (Morrell and Spielman, 1997; Ward et al., 2002), several techniques to suppress local distortions in EPI have been reported. In general, these standard approaches require an additional reference scan, for example, anatomical scans for imagebased methods (Kybic et al., 2000; Studholme et al., 2000), field maps (Jezzard and Balaban, 1995) or point spread function (PSF) maps (Zaitsev et al., 2004; Zeng and Constable, 2002). Additional scan time demands range from seconds in EPI-based field map approaches (Hutton et al., 2002; Reber et al., 1998) to several minutes in multi-echo reference acquisitions (Chen and Wyrwicz, 2001) or PSF mapping (Zeng and Constable, 2002) procedures. Moreover, computational time for the distortion correction can exceed several minutes for field mapping procedures or at least half an hour for some image-based approaches (Andersson et al., 2001; Studholme et al., 2000). The increase in acquisition and computational time makes it difficult or impossible to use the techniques for real-time applications. Here, we present a technique which corrects distortions along the PE direction in a single shot. It combines multi-echo EPI data acquisition with alternating PE gradients. This allows for simultaneous acquisition and distortion correction of the images. No additional reference scans are required and computation is efficient. We discuss the theory, implementation, and applications of this single-shot distortion correction. The reduction in geometric distortions is assessed in phantom MRI and human fMRI. Effects on SNR and BOLD sensitivity in fMRI are investigated.

Theory and procedure Theoretical background The presented undistortion technique combines the distortion correction method using PE gradients with alternating polarity (Bowtell et al., 1994; Chang and Fitzpatrick, 1992; Morgan et al., 2004) with a single-shot multi-echo EPI acquisition scheme (Mathiak et al., 2002). The multi-echo EPI sequence acquires three images at successive echo times with alternating polarity of the PE gradient blips (see PE gradient G PE in Fig. 1a). Hence, the second image

1069

will exhibit distortions of opposite direction but equal magnitude as compared to the first and third images (see original images in Fig. 1b; Bowtell et al., 1994; Chang and Fitzpatrick, 1992; Morgan et al., 2004). Therefore, we consider the first and second images only, and the results for the first image apply to the third image as well. Depending on the polarity of the PE gradient blip (G PE), an object at the position x with field offset DB will be projected to the positions x 1 and x 2 in the first and second images, respectively: x1 ¼ F1 ð xÞ ¼ x þ dx;

ð1Þ

x2 ¼ F2 ð xÞ ¼ x  dx

ð2Þ

with dx ¼

R

DBDt

ð3Þ

t0 þDt

t0

GPE ðtÞdt

and Dt being the echo spacing within the EPI readout echo train and t 0 being the beginning of the first PE gradient blip. If the correspondence of the positions x 1 and x 2 is known, it is straightforward to determine the local shift dx and the true position x of the object. Thus, calculation of this relationship for each x 1 and x 2 determines the undistortion transformations F 11 and F 21. Moreover, distortions modify image intensities. The true intensity q is estimated from the two image intensities q 1 and q 2 of the two images using the Jacobian of the transformation: qð xÞ ¼ q1 ðx1 Þ

dx1 dx2 ¼ q2 ðx2 Þ : dx dx

ð4Þ

Chang and Fitzpatrick (1992) have shown that F 1 and F 2 are monotonic functions, if the gradient in the field inhomogeneity is smaller than the applied effective gradient, that is, for distortions in the PE direction:  d ð DBÞ   Z t0 þDt .     ð5Þ GPE ðtÞdt  Dt  b dx t0 In the present study, the effective gradient amplitude in the PE direction was about 270 AT/m. Based on previous studies, maximal gradients in field inhomogeneities may be expected to be of a comparable size at 3 T (compare Deichmann et al., 2003) but restricted to rather small brain areas. Assuming monotonicity of F 1 and F 2, the correspondence between x 1 and x 2 can be found by integrating the signal intensities of the two images q 1, q 2 along the PE direction: Z

x1 ¼F1 ð xÞ

l

q1 ðF1 ðnÞÞdn ¼

Z

x2 ¼F2 ð xÞ

l

q2 ðF2 ðnÞÞdn:

ð6Þ

This derivation holds true for images acquired with the same parameters except for the polarity of the PE gradient. For instance, it could be applied to sequential acquisitions of EPI (Bowtell et al., 1994; Morgan et al., 2004). In multi-echo EPI, images are acquired in a single-shot at different echo times and exhibit different contrasts. Therefore, the signal intensity of the images should be adjusted. To that end, a virtual image is constructed with the same distortions as the first and third images and contrast properties similar to the second image acquired at TE2. Assuming a mono-exponential decay of the MR signal across the three images, the intensity I¯2 of the virtual image can be

1070

N. Weiskopf et al. / NeuroImage 24 (2005) 1068–1079

N. Weiskopf et al. / NeuroImage 24 (2005) 1068–1079

determined from the intensities I 1 and I 3 of the first and third image at echo times TE1 and TE3: P ln ðI 2 Þ ¼



   TE3  TE2 TE2  TE1 d lnðI1 Þ þ d lnðI3 Þ: TE3  TE1 TE3  TE1

ð7Þ

We will refer to I¯ 2 as the weighted geometric mean of I 1 and I 3. Implementation and procedure The formal derivation was based on continuous functions and must be discretized for implementation of the distortion correction, since MR data are acquired digitally with a limited sampling rate. Intermediate values are linearly interpolated from the data both for the signal integrals (Eq. (6)) used to establish the correspondence between the images and for the reconstruction of the undistorted images from the acquired images (Eq. (4)). The following undistortion procedure is applied to each volume of the time series separately. It consists of six steps (see Fig. 1b for a flow chart): (1) Data acquisition. At least three images are acquired in a single-shot using multi-echo EPI. The images are encoded with alternating PE gradients and thus exhibit opposite distortions. (2) Contrast adjustment and virtual second image. A virtual second image (I¯ 2) is determined by the weighted geometric mean of the first and third images (Eq. (7); the intensities are incremented by one, in order to avoid taking logarithms of values of zero). This virtual image possesses a similar contrast as the acquired second image (I 2) but shows opposite distortions. (3) Spatial smoothing. The virtual and acquired second images are smoothed in-plane using a Gaussian kernel (FWHM = 6 voxels) to improve SNR. (4) Estimation of distortion transformation. For each PE line, the correspondence between the smoothed version of the virtual and the acquired second image is determined by comparison of piecewise integrated signal intensities (Eq. (6)). In order to avoid mismatched integrals, the integrals across the image are normalized to be equal. The undistortion displacement dx is determined as half of the opposite shift between the voxels of the two images, that is, dx = (x 1  x 2)/2 (Eqs. (1) and (2)). A displacement map is calculated for the coordinate grid of the odd (first and third) and the even (second) images. (5) Distortion correction of images. For each image separately, an undistorted image is calculated by undoing the artifactual shift as determined in step 4 and the artifactual intensity modification as determined by the Jacobian of the displacement transformation (Eq. (4)). The Jacobian is determined by the numerical derivative of the displacement. (6) CNR optimization. A weighted average of the three undistorted images is composed for fMRI time series analysis.

1071

Images are weighted according to the BOLD sensitivity curve with the function f w(TE) = (TEd e TE/T2* )2, optimizing CNR in regions of dropouts due to decreased weighting of the third image at the long echo time (Mathiak et al., 2002). A value of T 2* = 35 ms at 3 T was chosen for the apparent relaxation time to account for brain regions with a short T 2*, such as parts of the cerebellum (Erb et al., 2003). Alternative weightings can be applied (Posse et al., 1999), such as constant weighting or weighting with the function f w(TE) = TEd e TE/T2*. For offline analysis, moreover, an adaptive weighting can be applied by estimating local T2* values (Mathiak et al., 2004).

Materials and methods An online version of the presented distortion correction and CNR optimization was implemented on a 3-T whole-body scanner (Magnetom Trio, Siemens Medical, Erlangen, Germany). All experiments were performed with a volume head coil using a single-shot multi-echo EPI sequence with alternating PE gradients (see below for imaging parameters). Real-time data export allowed for access of the data by auxiliary tools such as TurboBrainVoyager (Brain Innovation, Maastricht, The Netherlands) which performs online fMRI analysis (Goebel, 2001; Weiskopf et al., 2004). The processing time—including all image processing steps like Fourier transform, distortion correction, and data export—was less than 130 ms per slice (64  64 matrix) on the standard image reconstruction computer (Intel Xeon, 1.7 GHz, RAM 2 GBytes; Intel Corp., Santa Clara, CA, USA; implementation using the Siemens image calculation environment ICE VA23). For the sake of computational speed, the online output was restricted to the CNR-optimized images. However, all images were stored in the standard image data base in parallel for offline processing (for this evaluation study, processing was only performed offline). Offline algorithms were implemented in Matlab 5.3 (The MathWorks, Natick, MA, USA) and provided displacement maps and distorted, undistorted, and CNR-optimized images for a further detailed analysis. Experiment 1: geometry phantom To assess the amount of distortion reduction, a geometry phantom filled with doped water including plastic tubes was scanned using multi-echo EPI (12 transverse slices, slice thickness = 4 mm, gap between slices = 1 mm, flip angle a = 908, three images at effective echo times TE = 28, 46, 70 ms, bandwidth BW = 2.4 kHz/pixel, bandwidth in PE direction BWPE = 42 Hz/pixel, PE direction left–right, field of view FOV = 225  300 mm2, matrix size 48  64, partial Fourier PE factor 7/8, acquired PE lines = 42). The second-order shim (xy) was voluntarily detuned yielding displacements up to several voxels. A standard gradient echo image (FLASH, 12 transverse slices, slice thickness = 4 mm, gap between slices = 1 mm, TR = 533 ms, a = 608, TE = 5.19 ms, BW = 260 Hz/

Fig. 1. Single-shot multi-echo EPI acquisition scheme and distortion correction. (a) Three images were acquired after a single radio frequency (RF) excitation pulse at the echo times TE1 = 28 ms, TE2 = 46 ms, and TE3 = 70 ms. Polarity of the phase-encoding (PE) gradient G PE was inverted for the second image as can be seen as an inverted amplitude of the phase blips. Slice selective excitation using a gradient along the slice direction ( G S) and signal readout with a frequency encoding gradient G RO were identical to conventional EPI. (b) Displacements due to field inhomogeneities were reversed for the inverted polarity of the PE gradient (first row of images). Adjustment of contrasts of the images by weighted geometric mean (second row of images) and mapping of the images by the integration method (visualized on an equidistant grid; third row of images) yielded undistorted images (fourth row of images). The undistorted images were averaged to increase contrast-to-noise ratio (CNR), that is, sensitivity to the blood oxygen level-dependent (BOLD) effect (fifth row of images).

1072

N. Weiskopf et al. / NeuroImage 24 (2005) 1068–1079

pixel, PE direction left–right, FOV = 225  300 mm2, matrix size 96  128) which was minimally affected by distortions was obtained for comparison of distorted and undistorted EPI with the phantom geometry. Experiment 2: motor activation Five right-handed healthy subjects (two females, aged 23–34 years) participated in the study to determine the improvement of similarity between functional EPI and structural images by the undistortion procedure and to assess effects on fMRI time series analysis. For the latter, subjects performed a finger opposition task with the right hand according to a block design paradigm (duration 35 s per block, six rest blocks interleaved with five activation blocks). Informed written consent was obtained from the subjects prior to the study in accordance with the local ethics committee. Subjects were scanned using single-shot multi-echo EPI (38 coronal slices, slice thickness = 4 mm, gap between slices = 1 mm, repetition time TR = 3.5 s, flip angle a = 908, three images at effective echo times TE = 28, 46, 70 ms, bandwidth BW = 2.4 kHz/pixel, bandwidth in PE direction BWPE = 42 Hz/pixel, PE direction left–right, partial Fourier PE factor 7/8, field of view FOV = 224  168 mm2, matrix size 64  48, acquired PE lines = 42). We chose the imaging volume to cover the entire brain to allow for comparison with structural whole head scans. A T1weighted structural image (MPRAGE, 1 mm3 isotropic voxels, matrix size 256  224, partitions = 160, TR = 2300 ms, TE = 4.38 ms, inversion time IT = 1100 ms, a = 88) was obtained for the superposition of functional activations and for the comparison of distorted and undistorted EPI with an anatomical reference. Distortion correction, preprocessing, and statistical analysis Displacement maps were determined for each repetition, and each of the three images was undistorted separately. Additionally, CNR-optimized images were calculated as the weighted average of the three undistorted images. Seven time series entered the further analysis: three time series consisting of the original images, three time series consisting of the undistorted images, and one time series consisting of the CNR-optimized images. These time series were analyzed separately using the identical statistical parametric mapping procedure (SPM99, Wellcome Dept. Imaging Neuroscience, London, UK) to allow for comparisons. The first four volumes of each time series were excluded from analysis to suppress T1 equilibration effects. Time series were corrected for motion (Friston et al., 1995). To account for drifts and autocorrelations of the fMRI signal, the time series of each voxel was highpass- (cutoff period 140 s) and low-pass-filtered (Gaussian kernel, FWHM = 4 s). No spatial smoothing was applied. A general linear model (GLM) was applied to the time series of each voxel (Worsley and Friston, 1995). Activation blocks were modeled as a boxcar reference function which was convolved with a canonical hemodynamic response function (bigamma function with 6-s peak delay; Friston et al., 1998) to account for neurovascular coupling. Activated voxels were determined by testing for significant positive correlation of the measured signal with the modeled response, that is, significant signal increase due to the motor task. Conservative statistics were based on t statistics for each voxel and Bonferroni-corrected for multiple comparisons. Activations passing a fixed voxelwise threshold of t N 5.29 were considered significant (approximating a Bonferroni-corrected P b 0.01).

Effect of distortion correction on similarity of EPI and anatomy To assess the improvement in similarity of EPI images and structural images, the original and the undistorted first image (at TE = 28 ms) were each nonlinearly coregistered with an bindividualizedQ EPI template which was in turn matched to the individual anatomy using the structural image without the artifactual distortions. Coregistration was performed using the normalization procedure provided by the SPM99 package (Ashburner and Friston, 1999). The nonlinear components of this coregistration transformation in the PE direction can be attributed to artifactual image distortions. For quantification, the nonlinear residuals of the displacements were squared, summed, and normalized across the bindividualizedQ template brain mask. For the group of subjects, distortions in the first image were tested for decrease due to undistortion ( P b 0.05, Wilcoxon signed rank, one-sided). Due to heating of the gradient system and ferromagnetic shim elements, the shim and global frequency in the magnet may change during an imaging experiment. In EPI imaging, global shifts in frequency yield an image shift along the PE direction which can be usually observed as a slow linear drift in the motion parameter for the PE direction. To assess the dynamic correction of this frequency offset, the linear drift in this motion parameter was tested for a decrease due to the dynamic undistortion procedure ( P b 0.05, Wilcoxon signed rank, one-sided). Effect of distortion correction and CNR optimization on BOLD sensitivity Effects of distortion correction and CNR optimization on BOLD sensitivity were estimated from two types of measures. First, activations due to the motor task were evaluated. The sizes of the largest cluster of significantly activated voxels in the primary motor cortex and the cerebellum and the t values of the maximally activated voxel in the primary motor cortex and the cerebellum were determined. Across all subjects, these values were tested for an increase due to the undistortion and CNR optimization procedure as compared to the original images at three different echo times ( P b 0.05, Wilcoxon signed rank, onesided). Second, a measure of temporal SNR assessed improvement in CNR independently from functional activations. The temporal SNR was estimated from the GLM applied to the functional data; that is, the average signal was divided by the residual standard deviation of the fitted GLM. We determined the distribution of temporal SNR across all brain voxels. For the group of subjects, we tested the median of this distribution for an increase due to the undistortion and CNR optimization procedure ( P b 0.05, Wilcoxon signed rank, one-sided).

Results Experiment 1: geometry phantom The arbitrarily detuned second-order shim (xy) resulted in prominent nonlinear distortions (Fig. 2b). The EPI differed from the FLASH image which served as geometric reference, and the direction of distortions was dependent on the PE gradient polarity. Distortion correction reduced the artifactual distortions as shown in

N. Weiskopf et al. / NeuroImage 24 (2005) 1068–1079

1073

Fig. 2. Single-shot multi-echo EPI images of geometry phantom with arbitrarily detuned xy shim. (a) A slice of the imaged volume is presented. The contours of the gradient echo (FLASH) image are superimposed in red as a geometric reference on all images. (b) The three images (acquired at echo times TE = 28, 46, 70 ms) are highly distorted as compared to the real geometry measured by FLASH imaging whose contours are plotted in red. The images exhibit opposite distortions, because the phase-encoding (PE) gradient was inverted for the second image. (c) After distortion correction of each image, the match of the images and the match with the real geometry improved. Note that distortions were also reduced in areas close to signal dropouts as can be seen in the lower part of the third image at TE = 70 ms.

Fig. 2c. In areas close to signal dropouts at TE = 70 ms, the procedure improved the similarity of the images as well. Experiment 2: motor activation Effect of distortion correction on similarity of EPI and anatomy Distortion correction reduced the deformations (with the lower performance in some areas suffering from extensive dropouts and low SNR) and increased the similarity of functional EPI images with the anatomical reference image (see MPRAGE in Fig. 3). Nonlinear geometric distortions decreased by 34–44% (Fig. 4) yielding a significant effect across all subjects (Wilcoxon Z = 2.02; N = 5; P b 0.05). The reduction of linear drifts over time along the PE direction was 66–95% (Wilcoxon Z = 2.02; N = 5; P b 0.05); that is, effects of global frequency shifts due to heating of the gradient system and ferromagnetic shim elements were diminished. Effect of distortion correction and CNR optimization on BOLD sensitivity In all time series and subjects, task-correlated signal increases were observed within the motor cortex and the cerebellum (Fig. 3). The extent of the largest cluster within the cortex increased after distortion correction and CNR optimization as compared to the three original images (for all three comparisons Wilcoxon Z = 2.02; N = 5; P b 0.05; see Fig. 5a). The extent of the largest cluster within the cerebellum increased as well (comparison with the first original image: Wilcoxon Z = 2.02; second image: Wilcoxon Z = 1.75; third image: Wilcoxon Z = 1.75; N = 5; all P b 0.05; see Fig. 5b). Peak t values within the cerebellum were larger for the undistorted CNR-optimized images as compared to the original images (first and third images: Wilcoxon Z N 2.02; P b 0.05;

second image: Wilcoxon Z = 1.21; P N 0.11, n.s.; N = 5; see Fig. 5d). Peak t values within the cortex did not differ significantly (first image: Wilcoxon Z = 0.67; P N 0.25; second image: Wilcoxon Z = 1.21; P N 0.11; third image: Wilcoxon Z = 0.13; P N 0.44; N = 5; see Fig. 5c). The temporal SNR increased after distortion correction and CNR optimization as compared to any of the original images, that is, the median of the SNR across the whole brain increased (for all three images: Wilcoxon Z N 2.02; N = 5; P b 0.05; see Fig. 6).

Discussion We have presented the theory and implementation of a singleshot distortion correction for multi-echo EPI. For fMRI, the dynamic distortion correction was combined with weighted averaging of the images at different echo times in order to optimize sensitivity to the blood oxygen level-dependent (BOLD) effect. Processing times of less than 130 ms per slice on standard hardware allowed for its application in real-time fMRI, for example, for online neurofeedback (Weiskopf et al., 2004). Distortions were reduced, and the BOLD sensitivity increased as compared to images at single echo times. Improved localization and sensitivity in fMRI The procedure effectively reduced geometric distortions in phantom imaging (Fig. 2) and in functional imaging of human subjects (Figs. 3 and 4). Even in the case of severely distorted EPI images in the phantom due to a detuned shim, geometric distortions were reduced substantially (Fig. 2). Quantitative

1074

N. Weiskopf et al. / NeuroImage 24 (2005) 1068–1079

Fig. 3. Coronal slice through the cerebellum presenting the anatomical reference scan (MPRAGE; first column), the original images (second and third column), and the undistorted images (fourth to last column). For all subjects, undistorted images exhibited fewer geometrical distortions and improved similarity with the structural image, especially in the cerebellar region. To illustrate the effects of distortion correction on fMRI results, statistical analysis was performed separately for each image type. The results of the analyses are superimposed as statistical maps over the image types they were based on (corrected P b 0.01; extent N three voxels). Localization of functional activations resembled the expected neuroanatomical regions better after distortion correction. Furthermore, weighted averaging of the undistorted images improved the functional sensitivity as can be seen in the more prominent clusters of activated voxels in the CNRoptimized images (last column) as compared to the original images (second and third columns). Note that, for clearer and more detailed display, the third image (acquired at TE = 70 ms) is omitted and that the color scale may vary.

analysis of in vivo data showed that undistorted EPI images matched on average better to the structural brain images and thus to the individual anatomy (Fig. 4). In fMRI, shifts in global frequency over time are often observed, since the extensive use of the gradient system for EPI results in heating of the ferromagnetic shim elements and changes in the shim. The dynamic distortion correction reduced the spatial drift of the images due to changes in global frequency by 66–95%. The undercompensation is due to

the conservative character of the algorithm; that is, distortions are underestimated with increased homogeneous noise in the images. As a further improvement, an a priori rigid-body compensation of this shift can be applied. The global off-resonance effects might also be compensated online by using phase information from a navigator echo (e.g., Pfeuffer et al., 2002). Any procedure using dynamic parameters can possibly decrease the signal-to-noise ratio (SNR) of the time series (e.g.,

N. Weiskopf et al. / NeuroImage 24 (2005) 1068–1079

1075

after distortion correction and contrast-to-noise ratio (CNR) optimization (Figs. 3 and 5). An overall increase in temporal SNR across the brain was observed (Fig. 6), supporting the results obtained for the localized functional activations. Moreover, as expected from the reduction of distortions, the activations matched better to the locations predicted by functional neuroanatomy; for example, the overlap of activations with gray matter was increased (Fig. 3). Technical considerations and limitations

Fig. 4. Nonlinear distortions before and after undistortion of the first original image (acquired at TE = 28 ms) averaged across the brain. The first image was selected for analysis, because it provided the highest signal-tonoise ratio (SNR) and was most similar to the template used for the analysis and warping procedure.

motion correction, Mathiak and Posse, 2001). As the most important finding for fMRI applications, we showed in a study of motor activation an increased sensitivity to the BOLD effect

The presented method has certain limitations which are mainly due to decay of the magnetic resonance (MR) signal, that is, transverse relaxation, and the assumption that the measured signal is independent from the PE gradient polarity: (1) Differences in echo times of the multiple images should be kept to a minimum, because the images are adjusted to the same contrast by a weighted geometric mean only (Eq. (7)). In a strict sense, this assumption is valid only for a mono-exponential signal decay and for a sufficient SNR in all images. Especially in regions affected by severe susceptibility artifacts, considerable deviations from the exponential signal decay might occur (Kru¨ger et al., 2001) and yield wrong contrast adjustments. The consequential mismatch in the signal might lead to misregistrations. As shown by Chang

Fig. 5. Analysis of functional activations during a finger opposition task with the right hand. Statistical parametric mapping was performed for the original images and after the undistortion and contrast-to-noise ratio (CNR) optimization. For all subjects, statistical maps were thresholded at the fixed t value N 5.29 approximating a corrected P b 0.01. The extent of the largest clusters of contiguous voxels passing the threshold increased significantly after CNR optimization and distortion correction for (a) the cortex and (b) the cerebellum. The t values of the peak activations observed in (c) cortex and (d) cerebellum were comparable for all images.

1076

N. Weiskopf et al. / NeuroImage 24 (2005) 1068–1079

Fig. 6. Changes in temporal signal-to-noise ratio (SNR) due to the distortion correction and optimization of contrast-to-noise ratio (CNR). The temporal SNR was determined for each voxel inside the brain, and the median of this distribution across the brain was calculated. Medians of the SNR distribution of the original images and the undistorted CNR-optimized images are plotted.

and Fitzpatrick (1992), this misregistration is propagated through the image, but its amplitude is decreased in areas with higher signal [see Eq. (15) in Chang and Fitzpatrick, 1992]. Since these deviations from the exponential decay usually cause fast signal decay and low signal, the displacement maps might be inaccurate close to the susceptibility artifacts, but would be more accurate in the rest of the brain showing higher signal. The impact of small regions exhibiting non-exponential signal decay is further reduced by the spatial smoothing of the images used for estimation of the distortions. Furthermore, deviations from mono-exponential signal decay occur in regions containing compartments with differing T2* values, such as brain parenchyma and cerebrospinal fluid (CSF). However, the error caused by this partial volume effect and multiexponential decay does not exceed the error due to image noise in EPI (with SNR of approximately 100:1). The number of PE lines can probably not exceed a limit of 64 in gradient echo EPI at 3 T. Otherwise, the SNR might be too small for the third echo after encoding too many k-space lines. At 1.5 T, more PE lines might be acquired due to the longer apparent transverse relaxation time T2*. Therefore, spatial resolution along the PE direction is limited. Indeed, we reduced the number of PE lines by application of partial Fourier acquisition and a rectangular field of view (FOV). Parallel imaging using phased array coils might overcome this limitation by reducing the number of PE lines required for higher resolutions. It might even offer the possibility of acquiring more than three images after one radio frequency (RF) excitation pulse, thus further improving the efficiency of the distortion correction. (2) The distortion correction is based on a modification of the method by Chang and Fitzpatrick (1992) which uses inverted PE gradients. It is assumed that the signal is conserved, if all imaging parameters are fixed except for the PE gradient polarity. However, signal dropouts due to susceptibility artifacts depend on the direction of the PE gradient in EPI and might lead to a signal mismatch (Deichmann et al., 2002; Kannengiesser et al., 1999). Thus, differential rephasing due to field gradients in the PE direction can cause spurious displacements. However, the signal dropout might be predicted from the spatial shifts and thus might be included in the distortion correction by an intensity correction of

the original images or constraints on the smoothness of the displacement maps. Moreover, if the edges of the imaged object are not correctly detected, integration is based on different parts of the object for the different images and a signal mismatch might occur. Our approach did not incorporate the detection of the edges of the object. We assume an adequate SNR and that the edges of the field of view are approximately aligned. Violations of this assumption might yield misregistrations propagating through the image but with decreased amplitudes in areas with higher signal (Chang and Fitzpatrick, 1992). Since the main signal contribution is from inside the object and usually exceeds noise levels several times, the displacement maps might be inaccurate outside the object but would be more accurate inside the object. In general, any signal contribution from outside the object should be minimized and should be the same for all images at different echo times. Therefore, noise levels outside the imaged object should be similar for all images to avoid bias, because noise in the magnitude images is rectified and does not cancel out with spatial smoothing (Rayleigh distributed noise; Kannengiesser et al., 1999). In principle, this could be avoided by using complex valued images instead or by taking into account phase information. However, this poses practical problems, because the signal integral is not monotonous for complex values and the phase suffers from aliasing (Kannengiesser et al., 1999). Consequently, ghosting in images should be minimized (Morgan et al., 2004), and SNR should be similar for all images. In order to reduce the spacing between the echo times of the images and to acquire images with little ghosting, the presented method requires a high-performance gradient system. In this study, maximum gradient power was used, providing the highest effective bandwidth not exceeding safety limits for potential nerve stimulations in a body gradient system (for this imaging protocol: maximal slew rate d max c 200 mT/m/ms maximal gradient amplitude g max c 25 mT/m). Interestingly, even at this bandwidth and short readout time, distortions might still extend over several voxels as can be seen in the cerebellum (Fig. 3). Conventional EPI acquisitions usually employ smaller bandwidths to increase SNR which would result in even larger geometric distortions.

N. Weiskopf et al. / NeuroImage 24 (2005) 1068–1079

Alternative processing strategies The distortion correction might be improved by alternative algorithms for the matching of the images which are acquired with opposite gradient polarity yielding opposite distortions. In particular, dynamic time warping (DTW; Kannengiesser et al., 1999), coregistration based on mutual information (MI; Reinsberg et al., 2003), and least-square fits of the distortions using discrete basis functions (Andersson et al., 2003) have been proposed as alternatives to the integration-based method (Bowtell et al., 1994; Chang and Fitzpatrick, 1992; Morgan et al., 2004). DTW uses the information of the integrated signal along a PE line as well, but in contrast, it allows for easy inclusion of additional constraints on the displacement map. For instance, corrections for asymmetric dropout effects which lead to mismatch of the integrals (see above) may be included more easily into this algorithm. Reinsberg et al. (2003) reported that an MI coregistration in combination with a thresholded Sobel filter for detection of the edges of the object performed better in regions of signal dropouts and air cavities when compared to DTW and the integration methods. Andersson et al. (2003) showed that displacements can be modeled using a set of discrete basis functions. The parameters could be simultaneously estimated in 3D in a least-square sense which should increase the robustness. The method increased geometric fidelity in diffusion tensor imaging (DTI). DTW and MI coregistration have only been used to correct distortions in multi-shot gradient echo images and not for EPI. All three approaches have not yet been used for fMRI. Thus, the effect on temporal SNR was not studied. For application of any dynamic correction in fMRI, its effects on temporal signal stability and noise should be assessed, since it might increase temporal noise which can differ significantly from instantaneous image noise. Moreover, the demands these algorithms place on computational resources preclude real-time applications with current workstations. In the present study, intrinsic comparisons of undistorted and original images showed reduction of distortions and increase of BOLD sensitivity. In further studies, the performance of this technique could be compared with other undistortion and advanced shimming techniques such as point spread function mapping (Zeng and Constable, 2002; Zaitsev et al., 2004), field map methods (Cusack et al., 2003; Jezzard and Balaban, 1995; Jezzard and Clare, 1999), and real-time autoshimming (Ward et al., 2002). Moreover, effects of systematic variations of the parameters on the correction could be investigated; for example, the spatial smoothing of the original images might be varied, and spatial smoothing of displacement maps could be included. Additionally, alternative interpolation methods for calculation of the undistorted images from the distorted images might be studied (Morgan et al., 2004). Alternative techniques for distortionless echo-planar imaging Continuously updated EPI-based field maps have been proposed for dynamic distortion correction (Hutton et al., 2002; Jezzard and Clare, 1999). However, effects on time series were studied only for the double-shot technique which acquires two successive phase images and reduces the effective time resolution (Hutton et al., 2002). Recently, an approach based on spiral-in/spiral-out EPI acquisition was proposed for single-shot distortion correction

1077

(Sutton et al., 2004), which employs an iterative estimation of the signal equation. However, the iterative procedure results in an increased computational load and may not be applicable in a realtime application. Priest et al. (2004) have presented an alternative imaging technique (TRAIL) which uses a stimulated echo sequence to create a spatially interleaved spin echo and stimulated echo image, reducing the effective readout time and distortion. In comparison with multi-echo EPI, TRAIL results in an increased specific absorption rate (SAR) in the subject and spin system saturation. Thus, SAR safety limits might be exceeded, and the SNR might be decreased. Specific considerations for fMRI This study assessed geometric distortions due to field inhomogeneities and off-resonance effects. In fMRI, however, local signal dropout due to field inhomogeneities affects regions such as the orbitofrontal cortex, amygdala, and hippocampal regions. This technique does not specifically compensate such dropouts, but the multi-echo EPI acquisition provides images at short echo times which suffer from fewer dropouts. In addition, the undistortion technique can be combined with such compensation techniques (Deichmann et al., 2002, 2003; Posse et al., 2003b). Posse et al. (2003b) have reported a compensation scheme which used multiecho EPI with additional compensation gradients. Thus, multi-echo EPI might be used to reduce dropouts and geometric distortions simultaneously. The acquisition protocol in this study was specifically optimized for fMRI. Echo times were chosen to be close to T2* of gray matter to improve BOLD sensitivity (Mathiak et al., 2002; Posse et al., 1999). Partial Fourier imaging and a rectangular field of view reduced the echo time spacing between images and allowed for close stacking of the first two images (TE1 = 28 ms, TE2 = 46 ms) around the optimal echo time as given by T2*—similarly to spiral-in/spiral-out acquisition techniques (Preston et al., 2004).

Conclusions Recently, MRI systems with higher field strengths (N1.5 T) have become more widespread, offering superior BOLD sensitivity in fMRI and an improved signal-to-noise ratio (SNR). In particular, the improved SNR and BOLD sensitivity allows for functional imaging in real time, giving rise to new possibilities for quality assurance and assessment of neuronal processes; for instance, specific functional areas of the brain can be delineated online (Goodyear et al., 1997; Weiskopf et al., 2004), the subject’s BOLD response can be directly validated, and online feedback of the BOLD signal facilitates self-regulation of brain activity (deCharms et al., 2004; Weiskopf et al., 2003, 2004). However, geometric distortions in echo-planar imaging (EPI) increase with field strength and need to be compensated for. Most distortion correction procedures require additional time-consuming reference scans and in some cases lengthy computational processing, making it complicated or impossible to include these techniques in real-time imaging. Moreover, due to the static nature of the reference scan, they cannot correct for distortions that vary over time. The presented technique provides simultaneous dynamic correction of geometric distortions and optimization of BOLD sensitivity as validated in an fMRI experiment. It does not require

1078

N. Weiskopf et al. / NeuroImage 24 (2005) 1068–1079

additional data acquisition or extensive computation, allowing for its use in real-time fMRI.

Acknowledgments We thank M. Erb, B. Kardatzki, and W. Grodd for technical support and helpful discussions. We would also like to thank M. Borutta for the excellent technical assistance, N. Zink for revising the manuscript, and R. Deichmann for helpful discussions. The support with sequence programming by Siemens Medical (Erlangen, Germany) and in particular by S. Thesen is gratefully acknowledged. This study was supported by grants of the Deutsche Forschungsgemeinschaft (SFB 550/B5+B1 and Th 812/1-1) and BMBF. N.W. was supported by the WIN-research-group of the Heidelberg Academy of Science and Humanities (Germany).

References Andersson, J.L.R., Hutton, C., Ashburner, J., Turner, R., Friston, K., 2001. Modeling geometric deformations in EPI time series. NeuroImage 13, 903 – 919. Andersson, J.L.R., Skare, S., Ashburner, J., 2003. How to correct susceptibility distortions in spin-echo echo-planar images: application to diffusion tensor imaging. NeuroImage 20, 870 – 888. Ashburner, J., Friston, K.J., 1999. Nonlinear spatial normalization using basis functions. Hum. Brain Mapp. 7, 254 – 266. Bowtell, R.W., McIntyre, D.J.O., Commandre, M.-J., Glover, P.M., Mansfield, P., 1994. Correction of geometric distortion in echo planar images. Proc. SMR, 2nd Annual Meeting, Society of Magnetic Resonance, Berkeley, CA, USA, p. 411. Chang, H., Fitzpatrick, J.M., 1992. A technique for accurate magnetic resonance imaging in the presence of field inhomogeneities. IEEE Trans. Med. Imag. 11, 319 – 329. Chen, N.K., Wyrwicz, A.M., 2001. Optimized distortion correction technique for echo planar imaging. Magn. Reson. Med. 45, 525 – 528. Cusack, R., Brett, M., Osswald, K., 2003. An evaluation of the use of magnetic field maps to undistort echo-planar images. NeuroImage 18, 127 – 142. deCharms, R.C., Christoff, K., Glover, G.H., Pauly, J.M., Whitfield, S., Gabrieli, J.D.E., 2004. Learned regulation of spatially localized brain activation using real-time fMRI. NeuroImage 21, 436 – 443. Deichmann, R., Josephs, O., Hutton, C., Corfield, D.R., Turner, R., 2002. Compensation of susceptibility-induced BOLD sensitivity losses in echo-planar fMRI imaging. NeuroImage 15, 120 – 135. Deichmann, R., Gottfried, J.A., Hutton, C., Turner, R., 2003. Optimized EPI for fMRI studies of the orbitofrontal cortex. NeuroImage 19, 430 – 441. Erb, M., Klose, U., Grodd, W., 2003. T2* mapping for TE optimization in fMRI at different field strength. Proc. ISMRM 11, 1722. Friston, K.J., Ashburner, J., Frith, C.D., Poline, J.-B., Heather, J.D., Williams, S.C., Frackowiack, R.S.J., 1995. Spatial registration and normalization of images. Hum. Brain Mapp. 2, 165 – 189. Friston, K.J., Fletcher, P., Josephs, O., Holmes, A., Rugg, M.D., Turner, R., 1998. Event-related fMRI: characterizing differential responses. NeuroImage 7, 30 – 40. Glover, G.H., Li, T.Q., Ress, D., 2000. Image-based method for retrospective correction of physiological motion effects in fMRI: RETROICOR. Magn. Reson. Med. 44, 162 – 167. Goebel, R., 2001. Cortex-based real-time fMRI. NeuroImage 13, S129. Goodyear, B.G., Gati, J.S., Menon, R.S., 1997. The functional scout image: immediate mapping of cortical function at 4 Tesla using receiver phase cycling. Magn. Reson. Med. 38, 183 – 186. Hutton, C., Bork, A., Josephs, O., Deichmann, R., Ashburner, J., Turner,

R., 2002. Image distortion correction in fMRI: a quantitative evaluation. NeuroImage 16, 217 – 240. Jezzard, P., Balaban, R.S., 1995. Correction for geometric distortion in echo planar images from B0 field variations. Magn. Reson. Med. 34, 65 – 73. Jezzard, P., Clare, S., 1999. Sources of distortion in functional MRI data. Hum. Brain Mapp. 8, 80 – 85. Kannengiesser, S.A.R., Wang, Y., Haacke, E.M., 1999. Geometric distortion correction in gradient-echo imaging by use of dynamic time warping. Magn. Reson. Med. 42, 585 – 590. Krasnow, B., Tamm, L., Greicius, M.D., Yang, T.T., Glover, G.H., Reiss, A.L., Menon, V., 2003. Comparison of fMRI activation at 3 and 1.5 T during perceptual, cognitive, and affective processing. NeuroImage 18, 813 – 826. Krqger, G., Kastrup, A., Glover, G.H., 2001. Neuroimaging at 1.5 T and 3.0 T: comparison of oxygenation-sensitive magnetic resonance imaging. Magn. Reson. Med. 45, 595 – 604. Kybic, J., Thevenaz, P., Nirkko, A., Unser, M., 2000. Unwarping of unidirectionally distorted EPI images. IEEE Transactions on Medical Imaging 19, 80 – 93. Mathiak, K., Posse, S., 2001. Evaluation of motion and realignment for functional magnetic resonance imaging in real time. Magn. Reson. Med. 45, 167 – 171. Mathiak, K., Rapp, A., Kircher, T.T.J., Grodd, W., Hertrich, I., Weiskopf, N., Lutzenberger, W., Ackermann, H., 2002. Mismatch responses to randomized gradient switching noise as reflected by fMRI and wholehead magnetoencephalography. Hum. Brain Mapp. 16, 190 – 195. Mathiak, K., Hertrich, I., Grodd, W., Ackermann, H., 2004. Discrimination of temporal information at the cerebellum: functional magnetic resonance imaging of non-verbal auditory memory. NeuroImage 21, 154 – 162. Morgan, P.S., Bowtell, R.W., McIntyre, D.J.O., Worthington, B.S., 2004. Correction of spatial distortion in EPI due to inhomogeneous static magnetic fields using the reversed gradient method. J. Magn. Reson. Imaging 19, 499 – 507. Morrell, G., Spielman, D., 1997. Dynamic shimming for multi-slice magnetic resonance imaging. Magn. Reson. Med. 38, 477 – 483. Munger, P., Crelier, G.R., Peters, T.M., Pike, G.B., 2000. An inverse problem approach to the correction of distortion in EPI images. IEEE Trans. Med. Imag. 19, 681 – 689. Pfeuffer, J., Van de Moortele, P.-F., Ugurbil, K., Hu, X., Glover, G.H., 2002. Correction of physiologically induced global off-resonance effects in dynamic echo-planar and spiral functional imaging. Magn. Reson. Med. 47, 344 – 353. Posse, S., Wiese, S., Gembris, D., Mathiak, K., Kessler, C., Grosse-Ruyken, M.L., Elghahwagi, B., Richards, T., Dager, S.R., Kiselev, V.G., 1999. Enhancement of BOLD-contrast sensitivity by single-shot multi- echo functional MR imaging. Magn. Reson. Med. 42, 87 – 97. Posse, S., Binkofski, F., Schneider, F., Gembris, D., Frings, W., Habel, U., Salloum, J.B., Mathiak, K., Wiese, S., Kiselev, V., Graf, T., Elghahwagi, B., Grosse-Ruyken, M.L., Eickermann, T., 2001. A new approach to measure single-event related brain activity using real-time fMRI: feasibility of sensory, motor, and higher cognitive tasks. Hum. Brain Mapp. 12, 25 – 41. Posse, S., Fitzgerald, D., Gao, K., Habel, U., Rosenberg, D., Moore, G.J., Schneider, F., 2003a. Real-time fMRI of temporolimbic regions detects amygdala activation during single-trial self-induced sadness. NeuroImage 18, 760 – 768. Posse, S., Shen, Z., Kiselev, V., Kemna, L.J., 2003b. Single-shot T2* mapping with 3d compensation of local susceptibility gradients in multiple regions. NeuroImage 18, 390 – 400. Preston, A.R., Thomason, M.E., Ochsner, K.N., Cooper, J.C., Glover, G.H., 2004. Comparison of spiral-in/out and spiral-out BOLD fMRI at 1.5 and 3 T. NeuroImage 21, 291 – 301. Priest, A.N., Carmichael, D.W., DeVita, E., Ordidge, R.J., 2004. Method for spatially interleaving two images to halve EPI readout times: two reduced acquisitions interleaved (TRAIL). Magn. Reson. Med. 51, 1212 – 1222.

N. Weiskopf et al. / NeuroImage 24 (2005) 1068–1079 Reber, P.J., Wong, E.C., Buxton, R.B., Frank, L.R., 1998. Correction of off resonance-related distortion in echo-planar imaging using EPI-based field maps. Magn. Reson. Med. 39, 328 – 330. Reinsberg, S.A., Moore, L., Doran, S.J., Benekos, O., Leach, M.O., 2003. Mutual information-based correction of B0 inhomogeneity, susceptibility- and chemical-shift artefacts. Proc. ISMRM 11, 1049. Speck, O., Hennig, J., 1998. Functional imaging by I0- and T2*-parameter mapping using multi-image EPI. Magn. Reson. Med. 40, 243 – 248. Studholme, C., Constable, R.T., Duncan, J.S., 2000. Accurate alignment of functional EPI data to anatomical MRI using a physics-based distortion model. IEEE Trans. Med. Imag. 19, 1115 – 1127. Sutton, B.P., Noll, D.C., Fessler, J.A., 2004. Dynamic field map estimation using a single spiral-in/spiral-out acquisition. Magn. Reson. Med. 51, 1194 – 1204. Ward, H.A., Riederer, S.J., Jack Jr., C.R., 2002. Real-time autoshimming for echo planar timecourse imaging. Magn. Reson. Med. 48, 771 – 780. Weiskopf, N., Veit, R., Erb, M., Mathiak, K., Grodd, W., Goebel, R.,

1079

Birbaumer, N., 2003. Physiological self-regulation of regional brain activity using real-time functional magnetic resonance imaging (fMRI): methodology and exemplary data. NeuroImage 19, 577 – 586. Weiskopf, N., Mathiak, K., Bock, S.W., Scharnowski, F., Veit, R., Grodd, W., Goebel, R., Birbaumer, N., 2004. Principles of a brain–computer interface (BCI) based on real-time functional magnetic resonance imaging (fMRI). IEEE Trans. Biomed. Eng. 51, 966 – 970. Worsley, K.J., Friston, K.J., 1995. Analysis of fMRI time-series revisited—again. NeuroImage 2, 173 – 181. Yoo, S.S., Jolesz, F.A., 2002. Functional MRI for neurofeedback: feasibility study on a hand motor task. NeuroReport 13, 1377 – 1381. Zaitsev, M., Hennig, J., Speck, O., 2004. Point spread function mapping with parallel imaging techniques and high acceleration factors: fast, robust, and flexible method for echo-planar imaging distortion correction. Magn. Reson. Med. 52, 1156 – 1166. Zeng, H.R., Constable, R.T., 2002. Image distortion correction in EPI: comparison of field mapping with point spread function mapping. Magn. Reson. Med. 48, 137 – 146.