fMRI retinotopic mapping at 3 T: Benefits gained from correcting the spatial distortions due to static field inhomogeneity

Abstract

fMRI retinotopic mapping usually relies upon Fourier analysis of functional responses to periodic visual stimuli that encode eccentricity or polar angle in the visual field. Generally, phase estimations are assigned to a surface model of the cerebral cortex and borders between retinotopic areas are eventually determined following ad hoc phase analysis on the surface model. Assigning functional responses to a surface model of the cortex is particularly sensitive to geometric distortions of the 3D functional data due to static field inhomogeneity. Here, we assess and document the benefits gained from correcting the fMRI data for these effects, under standard experimental conditions (echo-planar imaging, 3.0-T field strength) and with well-chosen acquisition parameters (regarding slice orientation and phase-encoding direction). While it appears that, in the absence of correction, errors in the estimates of the borders between low-order visual areas do not then significantly exceed the variance of statistical origin, about half of the functional responses in a retinotopic experiment are misassigned to neighboring functional areas. Therefore, correction of the effects due to geometric distortions is important in any retinotopic mapping experiment and by extension in any fMRI experiment on the visual system.

Introduction

The human visual cortex is composed of a large number of functional areas with distinct functional properties (Felleman & Van Essen, 1991; Tootell, Tsao, & Vanduffel, 2003; Wandell, Brewer, & Dougherty, 2005). It is also organized hierarchically, to some extent. 

While V1, the lowest order visual area, is known to be centered within the calcarine sulcus and to extend over part of the banks from the surrounding gyri, the detection of anatomical landmarks such as the stripes of Gennari that would permit to delineate V1 requires applying a dedicated high-resolution imaging MR sequence, thus extending acquisition time typically by 45 min (Barbier et al., 2002). 

Beyond V1, there are no clear anatomical landmarks to precisely delineate any one among the visual areas. This delineation may be required in functional imaging studies, where the ability to localize functional imaging data with respect to the underlying functional architecture of the visual system may be important. 

With functional Magnetic Resonance Imaging (fMRI), it is possible to delineate a large number among the visual areas on the basis of their functional specificities. For establishing the boundaries of the lowest level visual areas, one relies upon their retinotopic organization—neurons from these areas respond to limited receptive fields, the centers of which are organized to form a continuous mapping between cortical surface and visual field. Importantly, at the border between adjacent, low-order, visual areas, the local representation of the visual field on the cortical surface changes its orientation—it is either a direct image or a mirror image of the visual field. The orientation of the representation of the visual field on the cortical surface is expressed as “the visual field sign”—which thus alternates between adjacent areas. These particular features of the low-order visual areas may be utilized to map retinotopy with fMRI. 

Stimuli for fMRI retinotopic mapping are designed to encode each position in the visual field by a unique pattern of temporal activation. While different approaches exist (Kraft et al., 2005; Slotnick & Yantis, 2003; Vanni, Henriksson, & James, 2005), fMRI retinotopic mapping most commonly relies upon the use of periodic stimuli that encode either eccentricity or polar angle in the visual field (Engel et al., 1994). 

To map eccentricity, slowly expanding or contracting annuli centered on the fixation point are used. To map polar angle, wedges that are positioned with the apex about the fixation point and slowly rotating in either clockwise or counterclockwise directions are used. These stimuli link each position along a visual field coordinate—eccentricity or polar angle—to a unique delay of the periodic stimulation of the neurons with the corresponding receptive fields (DeYoe et al., 1996; Engel, Glover, & Wandell, 1997; Hadjikhani, Liu, Cavanagh, & Tootell, 1998; Sereno, McDonald, & Allman, 1994; Wade, Brewer, Rieger, & Wandell, 2002; Warnking et al., 2002). As the BOLD response detected in fMRI represents a convolution of the neuronal response with the hemodynamic response function (Logothetis & Pfeuffer, 2004; Raichle, 1998; Ugurbil, Toth, & Kim, 2003), the delay measured represents the sum of the delay of interest and the hemodynamic delay. To cancel out the effects of the hemodynamic delay, fMRI responses are compared to stimuli that cover the visual field in opposite directions—i.e., for the polar angle visual coordinate, the neuronal delay is determined using wedges that rotate clockwise in a first run and counterclockwise in a second one, while for the eccentricity visual coordinate, the neuronal delay is obtained using annuli that expand in one run and that contract in another (Sereno et al., 1995; Warnking et al., 2002). Temporal delays are conveniently expressed in terms of phases in frequency domain, following Fourier transformation of the temporal signals. 

Data processing in fMRI retinotopic mapping thus differs from conventional fMRI in at least two respects. First, with a periodic stimulation paradigm and ensuing Fourier analysis, the main information of interest in fMRI retinotopic mapping is the phase of the signal in frequency domain (or, alternatively, its delay in temporal domain; amplitude is used only as a mask or weighting factor)—in conventional fMRI experiments, one is merely interested in the signal's amplitude; the response delay is supposed to be known. Second, to determine on the highly folded cortical surface where the borders between visual areas are located—i.e., where exactly on the cortical surface occurs the alternation of the visual field sign—signal phases need to be analyzed in their spatial context. While such analysis has been proposed in the 3D Cartesian space in which the data are acquired (Dumoulin et al., 2003), one usually proceeds to surface analysis of the functional responses. 

Surface analysis requires that an explicit model of the cortical surface be established first. Construction of the cortical surface relies upon the acquisition of high contrast 3D MR anatomical images, segmentation of the central layer of cortical gray matter on these images, removal of residual topological defects if any (Fischl, Liu, & Dale, 2001; Han, Xu, Braga-Neto, & Prince, 2002; Kriegeskorte & Goebel, 2001), and finally construction of a triangulated model of the cortex. For visualization purposes, this model may further be “inflated”(Fischl, Sereno, & Dale, 1999; Teo, Sapiro, & Wandell, 1997) or flattened (Grossmann, Kiryati, & Kimmel, 2002; Guérin-Dugué et al., 2000). 

In a second step, the fMRI responses to the visual stimuli need to be assigned to the surface model of the visual cortex. A difficulty arises because of the reduction in dimensionality—assigning volume-based data to the surface model is a non-homeomorphic process. Minor misalignment between anatomical and functional data may then induce very large errors in the surface representation of the functional data. As an example, misassignment of an activation to the opposite bank of the calcarine sulcus results in misassignment of the functional response to the opposite (upper vs. lower) visual hemi-field. The corresponding geodesic error along the cortical surface may then extend over as much as several centimeters, a distance over which the corresponding polar angle in the visual field may vary by several tens of degrees. Such variation also holds for the phase of the functional response measured in the retinotopic experiment. The issue of proper alignment of the functional and anatomic data needs to be addressed at the levels of data acquisition and of data preprocessing, prior to assigning the functional data to the surface model. A major cause for improper alignment of the functional and anatomic data sets—besides inter-scan subject motion—lies in the geometric distortion of the functional images due to static magnetic field inhomogeneity and low “apparent bandwidth” in the phase-encoding direction of the fast MR sequences usually applied in fMRI (single-shot echo-planar imaging, EPI; Jezzard & Balaban, 1995). Under typical conditions, this spatial distortion may amount up to several voxels. 

In the final step of the surface analysis, the borders of the low-order visual areas are determined from the detection of the changes of the visual field sign between adjacent areas (Sereno et al., 1995). To determine the visual field sign, it is then convenient to calculate the ratio of an oriented area measured using the local representation of the visual field coordinates with respect to some area measured using a locally isometric parametrization of the surface (this ratio is the Jacobian of the visual field representation on the surface). In a previous paper (Warnking et al., 2002), we referred to this quantity as the “visual field ratio” (VFR). The visual field sign is nothing else than the sign of the VFR. In addition, the visual area borders correspond to contour lines of zero VFR. Most of the spatial features of the VFR are present in the cortical representation of the polar angle coordinate in the visual field whose gradient reverses direction at the borders between visual areas. In contrast, the gradient of the cortical representation of the eccentricity coordinate in the visual field varies very smoothly across visual area borders. 

In an earlier paper (Warnking et al., 2002), we analyzed carefully all steps involved in fMRI retinotopic mapping, from stimulus design to the final steps of visual area delineation. We provided furthermore some suggestions for optimizing the stimuli used for mapping of polar angle retinotopy, for assigning the volume-based functional data to the cortical surface model, and for weighting phase information optimally to account for the local signal-to-noise ratio. The effect on the surface analysis of improper alignment of the functional and anatomic data, due to geometric distortions in the functional data, was not explicitly dealt with in that study, as experiments were performed at fairly low static field strength (1.5 T) and using an MR sequence (3D PRESTO) presenting reduced sensitivity to static field inhomogeneities—with this sequence, the differential distortion between anatomic and functional images is reduced by a factor of nearly five with respect to those expected with the single-shot EPI sequences traditionally applied in fMRI. Thus, the static field inhomogeneity was not a major issue in our earlier study. 

In the present study, we want to analyze and to document the benefits gained in fMRI retinotopic mapping—and more generally in any fMRI study on the visual system—from correcting the functional MR images for the geometric distortions, under standard experimental conditions: 3.0 T (or higher field) strengths, single-shot EPI MR acquisition of the BOLD contrasts, and surface analysis of the fMRI retinotopic data. 

Methods

Subjects

Experiments were performed on four healthy volunteers (25–45 years, 2 males, normal visual acuity), following project approval by the Institutional Review Board and written consent from the subjects. 

MR data acquisition

Images were acquired on a Bruker 3T Medspec S300 system with a gradient system providing a maximum gradient strength of 40 mT/m and a maximum slew rate of 133 T/m/s. A transmit–receive quadrature birdcage head coil (Bruker) was used. Structural, functional (retinotopic), and static magnetic field (inhomogeneity) data were acquired in a single scanning session. Static field homogeneity was optimized prior to acquiring the structural data, using the routine first-order shimming procedure. This preparation phase was not repeated prior to functional data acquisitions so that all images were acquired under identical static magnetic field distribution. 

Structural data were acquired in the sagittal orientation by means of a T1-weighted 3D MP-RAGE sequence with the following parameters: TR/TE/TI: 12/4.6/900 ms, excitation pulse angle: 8°, acquisition matrix: 256 × 256 × 176 (AP, CC, LR), isotropic nominal resolution: 1 mm, readout direction: CC, number of averages: 1, total measurement time of 14 min 40 s. 

Functional data acquired during retinotopic stimulus presentation were obtained using a 2D, gradient-recalled echo (GRE), multi-slice, single-shot EPI MR sequence with the following parameters: TR/TE: 2000/30 ms, excitation pulse angle: 77°, acquisition matrix: 72 × 64 (AP, LR), isotropic nominal resolution: 3 mm, bandwidth in the readout (LR) direction: 123 kHz, “bandwidth” in the phase-encoding (AP) direction: 1.7 kHz, 30 adjacent contiguous slices, thickness: 3 mm. Axial slices were angulated about the left–right axis to be approximately parallel to the calcarine sulcus. One of the central slices of the volume was positioned to contain as much of the calcarine sulcus as possible. Acquisition time per functional run was 7 min 50 s, allowing the acquisition of 235 volumes. 

If the local resonance frequency offset due to static field inhomogeneity is Δν and if δν represents the voxel “apparent bandwidth” in the phase-encoding direction, the signal displacement in units of voxel size is Δν/δν. At 3 T, the Larmor frequency (f ₀ = γB ₀/2π) for protons is 126 MHz. Thus, the signal displacement in the phase-encoding direction per ppm field offset, at 3 T, is 126 Δτ voxels/s, where Δτ is the overall measurement time in the phase-encoding direction. Δτ depends on the particular MR sequence used. Under our experimental conditions, this leads to a displacement of 4.7 voxels per ppm static field offset. 

Images representing signal phase delays due to off-resonance (field inhomogeneity) effects were derived from two 3D GRE MR sequences differing in the echo time only (5.4 and 14.5 ms, respectively) and with the following common parameters: TR: 25 ms, acquisition matrix: 64 × 256 × 64 (LR, AP, CC), acquisition duration: 2 min 25 s, nominal resolution: 4 × 1 × 4 mm³. Slice orientation was set identical to the functional images and high resolution was chosen in the phase-encoded direction of functional images where spatial distortion occurs. 

Stimuli for retinotopic mapping

Eccentricity in the visual field was mapped by a ring slowly contracting or expanding about the fixation point. Speed of contraction or expansion varied linearly with eccentricity. Thus, the activation wave on the cortical surface traveled at approximately constant speed (under the assumption of an exponential cortical magnification factor). When the ring reached maximum eccentricity (set at 8.4°), it wrapped around to be replaced by a new one at minimum eccentricity (set at 0.2°), and vice versa. Polar angle in the visual field was mapped by two wedges separated by a phase lag of 180° and rotating slowly at constant speed about the fixation point. Period of both eccentricity and polar angle stimulation was 32 s. Rings and wedges consisted of a black and white radial checkerboard flickering at a frequency of 4 Hz. Aspect ratio of the checks was kept constant by scaling their height linearly with eccentricity. The rationale underlying our choice of stimulation parameters may be found in Warnking et al. (2002). Stimuli were back-projected on a translucent screen positioned at the rear of the magnet. Subjects viewed this screen at a distance of about 2.2 m via a mirror fixed on the head coil. Stimulus images were generated under Matlab (Mathwork) and their timing was controlled using Presentation software (Neurobehavioral Systems). Successive stimulus images were presented at a frequency of 4 Hz, inducing the perception of an almost smoothly varying stimulus. Stimuli were started concomitantly with dummy MR excitations about 10 s prior to effective MR data acquisition so as to enable immediate response detection. Four retinotopic functional scans were acquired, one for each of the directions of motion of the rings and wedges. 

Modeling the cortical surface

Using a distributed Markovian method (Scherrer, Forbes, Garbay, & Dojat, 2009), the structural volume was segmented in three tissue types—cerebrospinal fluid, white matter, and gray matter. After the segmentation, the interface between volumes labeled white matter and cortical gray matter was extended to represent approximately the center of the cortical surface. Finally, some manual editing was performed to correct for topological defects detected in the occipital lobes. 

Only the cortical surface situated in the occipital lobes was then modeled. This allowed us eventually to flatten the surface model without the need for further cuts, thus preserving connectivity among visual areas. The region modeled was defined manually from the segmented volume, for each of the two hemispheres. It was delimited by two perpendicular planar cuts, the first one being positioned approximately parallel to and just posterior to the parieto-occipital sulcus and the second one being positioned approximately parallel to and about 3 cm antero-ventral of the calcarine sulcus. Within the delimited region, a triangulated model of the approximate center of the cortical surface interface was then generated following an approach based on the marching cubes algorithm (Lorenson & Cline, 1987). Distance between adjacent nodes of the triangulation was typically 1 mm. 

Flattened representation of the cortical surface model

Detailed information about the algorithm generating a flattened representation of the surface model may be found in Guérin-Dugué et al. (2000) and Warnking et al. (2002). Suffice it to mention here that we flattened the cortical surface for the purpose of display only. The flattening process is homeomorphic over the vast majority of the surface area and induces some distortion of local angles and distances. The analyses on the cortical surface that follow have been based entirely on the folded, triangulated surface model and are thus unaffected by these distortions. For the sake of display, we encode the local curvature of the cortical surface with different gray levels—gyri thus appear hyper-intense, while sulci appear hypo-intense on the gray scale used with the flattened representation of the cortical surface model. 

Volume reconstruction of retinotopic data

Functional echo-planar images were “deghosted” (Hennel, 1999), converted from proprietary scanner image format to SPM/Analyze, and preprocessed using SPM2 software. Two sets of functional images were considered. 

In the first set (“Uncorrected” EPI data set), the inhomogeneity of the static magnetic field was ignored—it was not taken into account during image reconstruction. Reconstructed images were realigned with respect to the first image of the series using rigid body motion correction (Ashburner, Neelin, Collins, Evans, & Friston, 1997). 

In the second set (“Corrected” EPI data set), the inhomogeneity of the static magnetic field was taken into account during functional image postprocessing. Using the phase information in the GRE images acquired at different echo times, a magnetic field map was first calculated (Cusack & Papadakis, 2002; Hutton et al., 2002). This was done with the SPM2 Fieldmap Toolbox. The magnetic field map was further used to compute a voxel displacement map and then to correct the functional images for the geometric distortions and to realign them with respect to the first one of the series. The conjoint field correction and realignment procedure was realized using the SPM2 Unwarp toolbox. In a final step, “uncorrected” and “corrected” EPI data sets were aligned with the structural data set using the SPM2 rigid body coregistration procedure based on mutual information metric. 

Volume-based analysis of retinotopic data

Following realignment with the structural data, the “uncorrected” and “corrected” EPI data sets were analyzed voxelwise by complex valued Fourier transformation, to determine amplitude and phase of the signals at the stimulation frequency. The hemodynamic response delay was accounted for by appropriately combining phase maps from stimuli moving in opposite directions (see details in Warnking et al., 2002). Each pair of functional time-series mapping one visual field coordinate gave thus rise to one pair of parametric data volumes: the phase of the stimulation at the visual field position represented locally and the response amplitude, representing the inverse of the uncertainty of the stimulation phase. Thus, four data volumes (two per visual field coordinate) were derived from each the “uncorrected” and the “corrected” EPI data sets (eight parametric volumes in total). 

The results of the volume-based analysis of the (four) retinotopic experiments were then stored, for each visual coordinate, as a pair of parametric data volumes containing, respectively, phase and signal-to-noise ratio at the stimulation frequency. Such four data volumes (two per visual coordinate) were derived from the “uncorrected” as well as from the “corrected” EPI data sets. In the following, we continue to express local visual field position in terms of the local phases of stimulation during expanding rings and clockwise rotating wedge stimuli, rather than polar visual field coordinates to which these phases could be converted, in order to keep relations to combine data originating from different voxels linear. 

Assignment of the phase information to the cortical surface model

Functional data were assigned to the cortical surface model for two distinct purposes: to generate retinotopic maps on the cortical surface model and to assess the distortion-induced misassignments. Two different approaches were used to assign data to the surface model, depending on the purpose of the maps obtained. 

When the purpose was generating retinotopic maps, phases were assigned to the cortical surface model as described in Warnking et al. (2002). Briefly, the phase at each node of the surface model was estimated as a linear combination of the phases stored in a selection of voxels from the parametric data volume (Equation 3 in Warnking et al., 2002). Voxels from the parametric data volume that were retained in this selection had their center at most 3 mm away from the closest node of the surface model and presented a functional response at the stimulation frequency with a signal-to-noise ratio exceeding the value of two. Weighting factors applied in the linear combination were determined as the product of the squared signal-to-noise ratio of the response in the voxel considered and a Gaussian weighting factor taking into account the Euclidean distance between the node from the surface model and the node closest to the voxel considered. Standard deviation of this Gaussian filter was set at 1.5 mm. To limit distance calculations, nodes presenting a geodesic distance exceeding 3 times the standard deviation of the filter were assigned a weight zero. This procedure was performed separately for the “uncorrected” as well as the “corrected” EPI data sets. The two retinotopic maps then obtained were used to qualitatively assess the superposition of the borders of the retinotopic areas as determined automatically by the method described in Warnking et al. (2002). 

When the purpose was assessing distortion-induced misassignments, we wanted to avoid the phase interpolations applied in Warnking et al.'s method to ensure unambiguous calculations of spatial shifts. We used only the parametric data volumes derived from the “uncorrected” EPI data set. We first rejected those data points in the “uncorrected” parametric data volumes that presented a signal-to-noise ratio at the stimulation frequency below 3 (corresponding to p > 10⁻³). We then applied the SPM2 Unwarp Toolbox to determine the corrected 3D coordinates from the data points retained in the “uncorrected” parametric data volume. For each visual coordinate, we were thus left with two 3D parametric data sets—one data set exhibiting phase information on certain nodes of a Cartesian grid (uncorrected data set) and one data set providing the same phase information on the corresponding nodes of a distorted grid (corrected data set). For the sake of simplicity, we will further refer to these parametric data volumes in a generic way (i.e., without mentioning explicitly to which visual coordinate they refer): PI_UNC (Phase Image on Cartesian grid—uncorrected data set) and PI_C (Phase Image on distorted grid—corrected data set). We further disregarded those data points in the two parametric data volumes that presented a distance to the closest node of the triangulated surface model exceeding 3 mm. Finally, the remaining data points were simply assigned to the closest node of the surface model. 

Generating retinotopic maps on the cortical surface model

Following assignment of the phase information to the surface model, retinotopic visual areas were automatically labeled and delineated on that model following the approach described in detail in Warnking et al. (2002). Briefly, this involved performing the following successive steps: (1) Application of a Gaussian filter (standard deviations of 3.5 mm and 7.0 mm in the polar angle and eccentricity phase directions, respectively) to reduce high-frequency noise of the phases assigned to the surface model; (2) calculation of the visual field ratio (VFR) map from the phase gradients with respect to a local two-dimensional coordinate system on the surface; (3) detection of candidate visual areas on the VFR map. Candidate areas are identified as contiguous areas exceeding a certain threshold on the absolute VFR map (∣VFR∣ > 8 (deg/mm²)) and on the SNR of the smoothed eccentricity and polar angle phase maps (SNR > 15); (4) selection of V1 as the largest among these candidates with negative VFR; (5) identification of the borders among retinotopic visual areas as the contour lines of zero VFR. Delineation of the visual areas from the VFR map was entirely performed in the two-dimensional Cartesian space of the flattened surface representation. 

Euclidean vs. geodesic distances

Distortion-induced misassignments to the cortical surface model were assessed by calculating the distances between corresponding nodes on the surface model derived from PI_UNC and PI_C. Euclidean distances were calculated straightforwardly. Geodesic distances along the cortical surface model were calculated using Dijkstra's shortest path algorithm (Dijkstra, 1959). 

Misassignments to opposite banks of gyri or sulci

If voxels are misassigned to a node from the cortical surface model while staying on the same bank of the gyrus or sulcus as the correctly positioned node, one expects distances between misassigned nodes and corresponding correctly positioned ones to be little sensitive to whether they are estimated along the surface model (geodesic distance) or directly in 3D (Euclidean distance). In contrast, when voxels are misassigned to an opposite bank from a gyrus or sulcus, the geodesic distance between misassigned and corrected positions will, in most cases, exceed significantly the Euclidean one. We therefore used the ratio R of the geodesic and Euclidean distances to identify cross-bank misassignments derived from PI_UNC—if R exceeded a threshold value set at 1.3, we considered the voxel having likely been misassigned to the wrong bank of a gyrus or sulcus. In what follows, we label these voxels as “cross-bank misassignments.” 

The rationale behind our choice of the R threshold was simply based upon modeling the bottom of a sulcus as a half circle. A voxel that is positioned at the bottom and that is wrongly assigned a quarter sector above would present an R value of 1.1. A voxel that is positioned at one top of the half circle and that is wrongly assigned to the opposite top (i.e., same level on the opposite bank of the sulcus) presents an R value of π/2 (see Figure 1). 

Figure 1

 

The bottom of a sulcus is modeled as a half circle. The blue and green voxels are assumed having been misassigned to the position of the gray voxel. They thus show an R ratio of π/2 and π/2 2 , respectively (geodesic distances are the lengths of the arcs, while Euclidean distances are the lengths of the chords).

Calculation of phase shifts

We then compared the phases of the nodes on the cortical surface model corresponding to cross-bank misassignments with the phases that would be expected from the phases at the surrounding nodes of the surface model—while disregarding all nodes corresponding to cross-bank misassignments. Phases at the nodes corresponding to cross-bank misassignments were estimated using a phase interpolation procedure similar to the one described earlier (procedure taking into account signal-to-noise ratios in the original 3D volumes and geodesic distances between the node considered and all nodes that contained phase information, other than those corresponding to cross-bank misassignments). Phase comparison relied simply upon subtracting the measured from the interpolated phases and by generating phase difference histograms. This was done for the phase maps derived from PI_UNC as well as for those derived from PI_C, and for each retinotopic visual coordinate. 

Misassignments to low-order visual areas

Following assignment of the phase information from the PI_C data sets to the surface model, we disregarded information in the nodes corresponding to cross-bank misassignments and then proceeded to phase interpolation and eventual delineation of the lowest order visual areas V1, V2v, V2d, V3v, and V3d as described earlier. Data points retained in PI_UNC were then projected onto the surface model, and for each of the five visual areas, we counted the number of data points that, following correction, remained either projected within the same visual area or were projected into another visual area or newly projected into the visual area considered. 

Results

Figure 2 shows histograms (all subjects) of the Euclidean and geodesic distances between corresponding nodes on the cortical surface model (i.e., nodes resulting from the assignment to the cortical surface model of corresponding voxels in PI_UNC and PI_C) and of the ratio R between these distances. Given the specific relationship between voxel shift and field homogeneity under our experimental conditions (4.7 voxels/ppm field offset), a Euclidean distance of 1 mm corresponds to a frequency shift of 9 Hz (i.e., 0.07 ppm). The distribution of Euclidean distances shows an FWHM of grossly 4 mm, corresponding to a frequency dispersion of 36 Hz. 

Figure 2

 

Histograms (all subjects and both hemispheres) of the distances between uncorrected and corrected voxel positions. (Left) Three-dimensional Euclidean distances. (Middle) Two-dimensional geodesic distances. (Right) Ratio R between geodesic and Euclidean distances.

Figure 3a shows a typical histogram (Subject 1) of the ratio R. In order not to overload Figures 3b, 3c, and 3d, only voxels with a geodesic displacement exceeding 8 mm are represented in this histogram. Figures 3b and 3c represent the correction-induced displacements on the planar surface model. In Figure 3b, we only considered voxels for which R was below 1.3. Most of these voxels remain on the same bank of gyri and sulci following correction. In Figure 3c, voxel shifts are shown for the voxels for which R exceeds the value of 2. A vast majority among these voxels is shifted across gyri or sulci to opposite banks. This may further be verified in Figure 3d representing these voxel shifts on the 3D surface model. The “cross-bank misassignments” (R > 2) represented 32% of the total number of voxels. 

Figure 3

 

(a) Typical histogram (Subject 1, right hemisphere) of the ratio R between geodesic and Euclidean distances. Only data points with a ratio R below 5 are represented in this histogram. (b) Voxel displacements following distortion correction. They are shown on a flattened representation of the cortical surface model (occipital lobe). Voxels considered are characterized by R ≤ 1.3. (c) Voxel displacements following distortion correction. They are shown on a flattened representation of the cortical surface model (occipital lobe). Voxels considered are characterized by R ≥ 2. (d) Same data as under (c). Displacements are, however, shown here with respect to the 3D cortical surface model. Many displacements cross gyri or sulci.

The two phase maps displayed on both sides of the phase color code strip in Figure 4 represent the (planar) surface model of a piece of cortical sheet encompassing part of the calcarine fundus for one subject (Subject 4). 

Figure 4

 

Planar surface model of a part of the cortical sheet encompassing the calcarine fundus (Subject 4; right hemisphere). Phases projected on the surface model correspond to the cortical representation of the eccentricity in the visual field. Phases on the left of the phase color code strip were projected from the uncorrected data set, while phases on the right of this strip were projected from a corrected data set. Large, square dots correspond to voxels with R ≤ 1.3, while small, diamond dots correspond to voxels with R ≥ 2. The two rectangular strips delineated by the black horizontal lines were drawn for the purpose of discussion (see text). The vertical axis corresponds to the antero-posterior direction.

Phases projected on the surface model correspond to the cortical representation of eccentricity in the visual field, as typically obtained. Phases on the left were projected from PI_UNC while those on the right were projected from PI_C. The vertical and horizontal directions in the figure correspond to first approximation to the direction along the calcarine fundus (direction along which the eccentricity in the visual field is encoded—also the direction into which MR signal spatial phase encoding was applied during the gradient-recalled echo EPI MR acquisition) and to the direction perpendicular to it (along which the polar angle in the visual field is encoded), respectively. Phases from voxels presenting a ratio R between geodesic and Euclidean displacements below 1.3 are represented as large square dots, while phases from voxels presenting a ratio R above 2 are presented as small diamond dots. Many among the latter must correspond to voxels having undergone cross-bank misassignment. This figure illustrates phase heterogeneity due to cross-bank misassignments in regions expected to present homogenous phase properties and the gain in phase homogeneity obtained following correction for the effects of static field inhomogeneity. In the example shown Figure 4 (Subject 4, right hemisphere), this is particularly striking in the cortical area delimited by the lower rectangle on the left sheet: well-positioned phases (R < 1.3; large square dots) correspond to values about 200° (yellow) while poorly positioned phases (R ≥ 2; small diamond dots) correspond to values about 175° (light green). Following correction for the effects of static field inhomogeneity (right sheet), most voxels with phase about 175° that were originally misassigned are now positioned at the level along the cortical sheet corresponding to the well-projected voxels with similar phase (upper rectangular area). 

Similarly, in Figure 5 (Subject 4, left hemisphere), following correction for geometric distortions (right sheet), more phases are projected and correctly positioned in the range between 100° and 250° and about 75°. 

Figure 5

 

Planar surface model of a part of the cortical sheet encompassing the calcarine fundus (Subject 4; left hemisphere). Phases projected on the surface model correspond to the cortical representation of the eccentricity in the visual field. Phases on the left of the phase color code strip were projected from the uncorrected data set, while phases on the right of this strip were projected from a corrected data set. Large, square dots correspond to voxels with R ≤ 1.3, while small, diamond dots correspond to voxels with R ≥ 2. The rectangular strips delineated by the black horizontal lines were drawn for the purpose of discussion (see text). The vertical direction corresponds to the direction parallel to the calcarine sulcus.

Figure 6 displays, for another subject (Subject 1), the surface models of the two hemispheres where phases projected correspond to the cortical representation of eccentricity in the visual field before and after geometric distortion correction. More voxels are correctly assigned to the 100° to 200° range on the surface model after spatial distortion correction. 

Figure 6

 

Planar surface model of a part of the cortical sheet encompassing the calcarine fundus for each hemisphere (left of the phase color strip: left hemisphere; right of the phase color strip: right hemisphere) for one subject (Subject 1). Phases projected on the surface model correspond to the cortical representation of the eccentricity in the visual field. Phases on the left of the vertical gray line were projected from the uncorrected data set, while phases on the right from this line were projected from a corrected data set. Large, square dots correspond to voxels with R ≤ 1.3, while small, diamond dots correspond to voxels with R ≥ 2. The rectangular strips delineated by the black horizontal lines were drawn for the purpose of discussion (see text). The vertical direction corresponds to the direction approximately parallel to the calcarine sulcus.

With our choice of phase-encoding direction (about parallel with the calcarine sulcus), geometric distortions due to field inhomogeneities are mainly in this direction into which eccentricity is encoded on the cortex. Thus, correction for these distortions will have little significant effect on the polar angle image (see Figure 7). 

Figure 7

 

Planar surface model of a part of the cortical sheet encompassing the calcarine fundus, left hemisphere, of two subjects (left: Subject 1; right: Subject 3). Phases projected on the surface model correspond to the cortical representation of the polar angle in the visual field. Phases above the phase color code strip were projected from the uncorrected data set, while phases below this strip were projected from the corrected data set. Large, square dots correspond to voxels with R ≤ 1.3, while small, diamond dots correspond to voxels with R ≥ 2. The horizontal direction corresponds grossly to the direction approximately perpendicular to the calcarine sulcus.

Figure 8 shows for all subjects how the borders among low-order visual areas compare when derived from PI_UNC (black limits) and from PI_C (white limits). Overall, Euclidean distances between corresponding borders determined using either data set may amount up to 4 mm on the folded surface model. 

Figure 8

 

Typical examples (top to bottom: Subject 1 to Subject 4) of how the borders among low-order visual areas compare when derived from uncorrected (black borders) and from corrected (white borders) 3D data sets. Borders were automatically delineated by the tool using the hypotheses that V1 is the largest area and that the visual field sign changes between adjacent areas. Manual intervention was sometimes required to refine borders, especially of V3.

Figure 9 represents, pooled for both hemispheres, the phase difference histograms averaged over all subjects both for eccentricity- and polar-angle-encoding stimuli. 

Figure 9

 

Phase difference histograms (all subjects). (Left) Eccentricity-encoded stimulus. (Right) Polar-angle-encoded stimulus. Red: derived from uncorrected data set. Green: derived from corrected data set. We considered reference points as points presenting no cross-bank misassignment, i.e., showing a low R (R < 1.3). Considering that the phase varies smoothly along the cortical surface, we used the phases of these reference points to determine by interpolation the phase for points presenting cross-bank misassignments, i.e., for which R ≥ 2. For these points, we then computed the difference between measured and interpolated phases.

In Figure 10, the phase difference histograms are detailed for each visual area for eccentricity-encoding stimuli. For recall, the phase differences taken into account are those from voxels presenting surface model misassignments characterized by R ≥ 2, thus including those from most cross-bank misassignments. Phase differences are those between measured phases and phases derived by interpolation from well-positioned data points (R < 1.3). Histograms shown in red were derived from PI_UNC while those shown in green were derived from PI_C. Figures 9 and 10 show that, after correction, for geometric distortions the phases assigned to “cross-bank voxels” have become less biased with respect to phases of neighborhood voxels (histograms centered to zero). 

Figure 10

 

Phase difference histograms (all subjects) for eccentricity-encoded stimulus in all delineated retinotopic areas (d: dorsal, v: ventral). Red: derived from uncorrected data set. Green: derived from corrected data set. We considered reference points as points presenting no cross-bank misassignment, i.e., showing a low R (R < 1.3). Considering that the phase varies smoothly along the cortical surface, we used the phases of these reference points to determine by interpolation the phase for points presenting cross-bank misassignments, i.e., for which R ≥ 2. For these points, we then computed the difference between measured and interpolated phases. This was done for the uncorrected (red) and corrected (green) data sets.

Figure 11 illustrates the misassignments to the surface model in the absence of corrections for static field inhomogeneity and the extent to which these misassignments induce errors in the assignment of functional responses to particular visual areas. The typical results shown here are those from Subject 2. The area borders have been determined using PI_C. In order not to overload the figure, only data points are shown for which the geodesic displacement following correction exceeds 10 mm. The red dots present their location on the planar surface model prior to correction and the green dots their position following correction. Yellow arrows represent the displacements along the surface model. The figure illustrates that, in the absence of corrections for static field homogeneity, many data are wrongly assigned to particular visual areas. 

Figure 11

 

Typical voxel displacements (Subject 2) following distortion correction. They are shown on a flattened representation of the cortical surface model (occipital lobes, left and right hemispheres), with delineation of the low-order visual areas. Voxels considered are characterized by a geodesic displacement exceeding 10 mm. Red dots: location prior to correction. Green dots: location following correction. Yellow arrows: displacements along the surface model. Prior to correction, many data are misassigned to neighboring visual areas.

Figure 12 summarizes, integrated over all subjects but separately for the left and right hemispheres, how the corrections for static field inhomogeneity affect assignment of the functional data to the low-order visual areas up to V3, as delineated using PI_C. The column charts show, per area, the number of voxels for which the correction did not affect the initial area assignment, the number of voxels that the correction assigned newly to the area, and the number of voxels that following correction were assigned to another area. The corresponding values for left and right hemispheres are shown, respectively, in Tables 1 and 2. 

Figure 12

 

Per visual area and per hemisphere (all subjects), the number of voxels for which the correction did not affect the initial area assignment (black), the number of voxels that the correction assigned newly to the area (light gray), and the number of voxels that following correction were assigned to another area (gray).

Table 1

 

Percentages of voxels correctly assigned, no longer assigned, and newly assigned to each area following distortion correction for left hemisphere (all subjects).

	V1	V2v	V3v	V2d	V3d	Overall
Not affected (%)	57	43	44	41	37	48
No longer assigned (%)	10	25	27	13	24	18
Newly assigned (%)	33	31	29	46	39	34

Table 2

 

Percentages of voxels correctly assigned, no longer assigned, and newly assigned to each area following distortion correction for right hemisphere (all subjects).

	V1	V2v	V3v	V2d	V3d	Overall
Not affected (%)	54	44	47	39	34	46
No longer assigned (%)	7	24	26	22	33	19
Newly assigned (%)	39	32	27	39	33	35

Discussion

In the presence of static field inhomogeneities, transverse magnetization in an MRI experiment acquires locally a phase offset in the so-called “rotating frame” (Schlichter, 1978). This phase offset is proportional to the signed local static field offset and the time of measurement of the transverse magnetization. Conventional spatial decoding in MRI then leads to spatial mispositioning and/or blurring of the information of interest according to the particular MR acquisition and reconstruction techniques applied. With a Cartesian MR acquisition technique such as EPI, spatial mispositioning is induced mainly in the direction into which spatial phase encoding has been applied during MR acquisition, leading to geometric distortions in this direction that are proportional to the local static field offsets. 

In this study, we analyzed and documented the potential benefits gained in surface analysis of retinotopic fMRI data from correcting the fMR images for the geometric distortions due to static field inhomogeneities. We assumed experimental conditions that have become standard in functional MRI, i.e., 3.0-T static field strength and gradient-recalled echo EPI MR acquisition. 

The need for dealing appropriately with the geometrical distortions due to static field inhomogeneity stems mainly from the need for analyzing phase properties of the functional responses to retinotopic stimuli in their spatial context along the cortical surface. The distance between corresponding voxels assigned to the cortical surface model prior to and following correction for the static field inhomogeneities may differ considerably according to whether one considers the Euclidean distance in 3D space or the geodesic distance along the cortical surface—this is particularly true for voxels misassigned to the opposite bank of a gyrus or sulcus on the highly folded cerebral cortex. Thus, if no correction for the geometric distortions is applied, voxels (mis)assigned to the cortical surface model may exhibit significant phase errors in fMRI retinotopic mapping. 

The significant difference between Euclidean and geodesic distances between voxels assigned to the surface model prior to and after distortion correction has been illustrated in Figure 2 using the data pooled from our set of subjects: while distortions cause positional errors of at most 11 mm in 3D space, the latter extend up to 40 mm along the cortical surface. Over such a distance along the cortical surface, the phase of the responses to retinotopic stimuli may vary considerably. Phase discrepancies due to misassignments on the cortical surface are further illustrated in Figure 4: the phases of the retinotopic responses from neighboring voxels assigned to the cortical surface model may differ by as much as 70 degrees if no correction for static field inhomogeneities is applied. These incongruences are to a large extent relieved following correction for the distortions. 

The correction for static field inhomogeneities sharpens significantly the phase difference histograms (Figures 9 and 10). In other words, phase match between directly measured values and values interpolated from the surrounding data improves significantly upon correction for the spatial distortions due to static field inhomogeneity, a further illustration of the benefits gained from correcting for the geometric distortions. Phase measurements are affected by noise. The observation that the distribution is not widened for the uncorrected data reflects the fact that the standard deviation of the phase differences due to field homogeneities is masked by the standard deviation due to the uncertainties in the phase measurements. 

The phase correction is more efficient when applied to data sets corresponding to eccentricity-encoding stimuli (see Figures 4, 5, and 6) than to data sets corresponding to phase-encoding stimuli (see Figure 7). Indeed, distortions are in the phase-encoding direction, i.e., here in the antero-posterior direction, which is roughly the direction of varying eccentricity onto the cortical ribbon according to retinotopy property. 

It turns out (Figure 8) that the errors of the area borders derived from uncorrected data sets are comparable to the variance obtained when assessing reproducibility of the delineation procedure using different data sets acquired under identical experimental conditions and processed similarly (Warnking et al., 2002). The limited sensitivity of the delineation to the geometrical deformations in this study is to be ascribed to the ad hoc choice of the spatial phase-encoding direction during MR acquisition of the functional data. Orientation of the image planes—transverse—was such that the central slice of the image volume encompassed most of the calcarine fundus. Importantly, the spatial phase-encoding direction was then chosen antero-posterior. Thus, the geometrical deformations due to static field inhomogeneities occurred mainly in the antero-posterior direction. This direction is that into which eccentricity in the visual field is encoded, while the polar angle in the visual field is encoded into a direction that is oriented perpendicularly to it—this general feature holds to certain extent for all low-order visual areas considered in this study. As the borders between visual areas are determined on the basis of the VFR, as most of the spatial features of the VFR are present in the cortical representation of the polar angle visual coordinate, and given the particular conditions of MR acquisition that were ours in this study, it did not come out as completely unexpected to us that the borders between low-order visual areas present limited sensitivity to static field inhomogeneities. 

The latter result does by no means imply that one may neglect to correct the functional data sets for the effects of static field inhomogeneities. As illustrated in Figure 11 for one subject, without correction many voxels are wrongly assigned to a particular visual area. Our results (Figure 12 and Tables 1 and 2) indicate that about half of the functional responses are reassigned to a particular visual area after correction for the effects of static field inhomogeneities. Here, these reassignments have occurred for an fMRI retinotopic mapping type of experiment. Similar rates of reassignments are to be expected for any cognitive fMRI experiment on the visual system performed under similar experimental conditions. While most of these reassignments are presumably correct, as indicated by quantitative analysis of Figures 4–7, it may not be excluded that some errors are introduced during processing. 

It is further noteworthy that alignment of functional and anatomical data in the occipital lobe may be improved using coregistration algorithms that focus on the region of interest (here the occipital lobe; Yassa & Stark, 2009). 

The present study was performed on an MR system equipped with a single detection channel and first-order shimming only. More recent MR systems provide multi-channel detection and higher order shimming. The use of multi-channel detection in combination with ad hoc data acquisition schemes (Griswold et al., 2002; Pruessmann, Weiger, Scheidegger, & Boesiger, 1999) permits to reduce significantly the readout time of the MR signal and the geometric distortions due to field inhomogeneities as a consequence. The use of higher order shimming may further reduce the field inhomogeneities over the volume of interest. Our results indicate that, with fMRI of the visual system, full advantage needs to be taken of these features when available. In the same vein, given the bandwidth and the field of view applied in the readout direction during acquisition of the functional data, the readout gradient strength was in the present study somewhat less than half the maximum gradient strength allowed by the system. We could thus have gained more than a factor of two in the geometrical distortions of the functional volumes—at the expense of a loss in signal-to-noise ratio. The particular choice of acquisition parameters was considered to represent an acceptable compromise between spatial distortions and signal to noise at the start of the study. 

Fully correcting for field inhomogeneities requires confronting various problems. The solution we propose suffers from two limitations: (1) Noise: measurements of the field shift certainly contain noise that hamper the computation of a “real field shift” and (2) Evolution: the field shift can change over time depending on head motion or changing oxygenation in large draining veins, then the computation of a single static correction is certainly limited in accuracy. 

Conclusion

In conclusion, we have assessed and illustrated the benefits gained from correcting the data in fMRI experiments on the visual system for the effects of static field inhomogeneities. We took advantage of the retinotopic properties of the low-order visual areas that permit to precisely determine the borders between adjacent areas on a cortical surface model following correction for the static field inhomogeneities. While, under our carefully optimized experimental MR acquisition conditions, it appears that the borders between adjacent areas present limited sensitivity to the static field inhomogeneities, about half of the functional responses to a retinotopic experiment are misassigned to a neighboring functional area. It may be anticipated that similar misassignments would also occur in any fMRI experiment on the visual system and more generally on any other cognitive system. 

Acknowledgments

Flor Vasseur is a recipient of a PhD grant from the Ministère de l'Enseignement Supérieur et de la Recherche (France). 

We would especially like to thank the two reviewers for their time and suggestions that have greatly helped improve the quality of the final paper. 

Commercial relationships: none. 

Corresponding author: Dr. Michel Dojat. 

Email: [e-mail]. 

Address: BP 170, Grenoble 38042, France.