Laminar analysis of 7 T BOLD using an imposed spatial activation pattern in human V1

doi:10.1016/j.neuroimage.2010.05.005

Neuroimage. Author manuscript; available in PMC 2011 October 1.

Published in final edited form as:

Neuroimage. 2010 October 1; 52(4): 1334–1346.

Published online 2010 May 9. doi: 10.1016/j.neuroimage.2010.05.005

PMCID: PMC3130346

NIHMSID: NIHMS213967

Laminar analysis of 7 T BOLD using an imposed spatial activation pattern in human V1

Jonathan R. Polimeni,¹ Bruce Fischl,^1,² Douglas N. Greve,¹ and Lawrence L. Wald^1,³

¹Athinoula A. Martinos Center for Biomedical Imaging, Department of Radiology, Massachusetts General Hospital, Harvard Medical School, Charlestown, MA, USA

²Computer Science and Artificial Intelligence Laboratory (CSAIL), Massachusetts Institute of Technology, Cambridge, MA, USA

³Harvard-MIT Division of Health Sciences and Technology, Massachusetts Institute of Technology, Cambridge, MA, USA

Corresponding author: Jonathan R. Polimeni, Ph.D., A.A. Martinos Center, MGH, Bldg 149 Thirteenth St., Suite 2301, Charlestown, MA 02129, Phone: (617) 306-3125, Email: jonp/at/nmr.mgh.harvard.edu

The publisher's final edited version of this article is available at Neuroimage

See other articles in PMC that cite the published article.

Publisher's Disclaimer

Abstract

With sufficient image encoding, high-resolution fMRI studies are limited by the biological point-spread of the hemodynamic signal. The extent of this spread is determined by the local vascular distribution and by the spatial specificity of blood flow regulation, as well as by measurement parameters that (i) alter the relative sensitivity of the acquisition to activation-induced hemodynamic changes and (ii) determine the image contrast as a function of vessel size. In particular, large draining vessels on the cortical surface are a major contributor to both the BOLD signal change and to the spatial bias of the BOLD activation away from the site of neuronal activity. In this work, we introduce a laminar surface-based analysis method and study the relationship between spatial localization and activation strength as a function of laminar depth by acquiring 1 mm isotropic, single-shot EPI at 7 T and sampling the BOLD signal exclusively from the superficial, middle, or deep cortical laminae. We show that highly-accelerated EPI can limit image distortions to the point where a boundary-based registration algorithm accurately aligns the EPI data to the surface reconstruction. The spatial spread of the BOLD response tangential to the cortical surface was analyzed as a function of cortical depth using our surface-based analysis. Although sampling near the pial surface provided the highest signal strength, it also introduced the most spatial error. Thus, avoiding surface laminae improved spatial localization by about 40% at a cost of 36% in z-statistic, implying that optimal spatial resolution in functional imaging of the cortex can be achieved using anatomically-informed spatial sampling to avoid large pial vessels.

Introduction

Functional MRI at high magnetic field strengths such as 7 T provides several advantages over conventional field strengths, such as increased T₂* contrast changes and increased image signal-to-noise ratio (SNR). The resulting increase in available functional contrast-to-noise ratio (CNR) can be exploited by acquiring data with reduced voxel sizes to produce finer spatial sampling of the functional signal across the brain. Furthermore, recent advances in both gradient coil and receive array coil technology can limit the vulnerability of high-field imaging to both image distortions [de Zwart et al., 2002] and signal drop-outs [Merboldt et al., 2000] due to macroscopic susceptibility gradients. The use of smaller voxels also brings functional imaging into a thermal noise dominated regime by reducing the contribution of physiological noise [Triantafyllou et al., 2005]. Thus high-resolution functional studies receive the full benefit of high-field magnets and highly-parallel receive array coils since these hardware-based advances boost signal levels relative to thermal noise. The use of smaller voxels also increases the proportion of voxels containing a large percentage of gray matter. With the feasibility of reduced voxel sizes, biological factors—such as the spatial specificity of the hemodynamic response—become the dominant limitation in the accurate localization of functional activation with BOLD imaging.

At high field strengths, the signal generated with standard BOLD contrast is relatively insensitive to intra-vascular changes due to the short T₂* of blood, and hence most signal changes are detected in the extra-vascular compartments [Yacoub et al., 2003]. The standard gradient-echo EPI used for fMRI acquisitions is broadly sensitive to extra-vascular signal changes surrounding vessels of all sizes above roughly that of the capillary, including those surrounding (relatively) large draining veins lying along the pial surface [Boxerman et al., 1995]. These pial veins pool deoxygenated blood from extended territories and therefore potentially displace the observed signal changes from the site of neuronal activity. The distribution of vessel diameters and density, however, is not uniform throughout the cortex. There is a hierarchy of vessel morphology arranged in vascular layers parallel to the cytoarchitectonic layers of the cortical gray matter. The highest capillary density lies within a vascular layer coinciding with Layer IV of the neocortex, and each vascular layer is punctuated by intermittent diving venules and arterioles oriented radially and which are in turn, connected to larger feeding arteries and draining veins coursing along the pial surface [Duvernoy et al., 1981; Lauwers et al., 2008]. Typical pial veins range in diameter from 1 to 4 mm, while diving venules running radial to the cortical lamination range in diameter from 20 to 150 μm, and capillaries within the capillary bed range in diameter from 3 to 8 μm [Duvernoy et al., 1981].

The significant contribution of deoxygenated hemoglobin carried downstream from the true site of neuronal activation by large vessels has motivated the use of alternative methods to standard gradient-echo BOLD imaging to increase spatial specificity in studies that require high spatial resolution. For example, other aspects of the hemodynamic response, such as changes in blood flow or blood volume, have been suggested to provide higher spatial specificity, yet these techniques suffer from reduced sensitivity [Ogawa et al., 1993; Bandettini et al., 1994; Boxerman et al., 1995]. It has also been suggested that by employing a high temporal sampling rate the “initial dip” immediately preceding the standard “positive” BOLD response could be detected, which may provide higher spatial specificity since it may reflect rapid deoxygenation prior to the increases in blood flow and volume that give rise to the positive BOLD response [Ernst & Hennig, 1994; Yacoub et al., 2001b]. However the initial dip effect is significantly smaller than positive BOLD response and can be inconsistent, as has been seen in monkey visual cortex [Logothetis, 2000]. In addition to acquisition strategies, post-processing techniques have also been employed, including the manual identification of larger veins and subsequent exclusion of voxels near to those vessels to reduce the number of contaminants that lower overall specificity [Yacoub et al., 2001a; Koopmans et al., 2009].

Spin-echo-based sequences have been shown to be preferentially sensitive to extra-vascular BOLD signal changes near smaller vessels, peaking at approximately the diameter of capillaries or small venules [Boxerman et al., 1995]. They are therefore advocated for high-resolution fMRI, particularly for high field strength studies where the sensitivity gap between spin-echo contrast and gradient-echo contrast experiments narrows [e.g., Duong et al., 2003; Yacoub et al., 2007]. Although spin-echo EPI techniques yield increased spatial accuracy, there are several drawbacks to this approach, such as the reduced size of the BOLD effect [Ogawa et al., 1993], increased Radio Frequency (RF) power deposition, and the sensitivity of the spin-echo signal strength to excitation flip angle inaccuracies arising from the shorter RF wavelength at high field.

Previous studies examining changes in fMRI activation as a function of laminar depth, both in humans [Gati & Menon, 2002; Ress et al., 2007; Koopmans et al., 2009; Yacoub et al., 2009] and animal models [e.g., Silva & Koretsky, 2002; Goense & Logothetis, 2006; Harel et al., 2006; Zhao et al., 2006; Moon et al., 2007; Smirnakis et al., 2007; Jin & Kim, 2008], have detected either earliest or greatest responses in the central cortical layers when all cortical layers along a cortical column are activated. This elevated Layer IV BOLD response could also arise from the increased blood volume in this lamina, without requiring laminar-level regulation of blood flow or oxygen consumption.

Existing approaches rely on analyses restricted to a small region of interest where the cortex is locally flat, and have resorted to hand-drawn profiles through cortical depths orthogonal to the gray matter boundaries in small, flat patches of cortex away from the cortical folds [Ress et al., 2007; Koopmans et al., 2009]. They are therefore not able to study patterns of activity at a particular cortical depth over larger areas of folded cortex. Here we introduce a new surface-based laminar analysis method that enables sampling of the fMRI signals within several intermediate cortical layers, and use it to quantify relative spatial specificity as a function of cortical depth or proximity to the draining veins. This surface-based laminar analysis is capable of analyzing patterns of activity within the cortex at a particular cortical depth over the full extent of the cortex, and is therefore not limited to a small region of interest or restricted to a flat patch of cortex exhibiting little or no curvature.

In this study, we sought to exploit the laminar organization of the venous anatomy by acquiring gradient-echo EPI data with small isotropic voxels, and limiting the analysis to voxels away from the large pial vessels within the central and deep cortical laminae, where the vessels are considerably smaller and more densely packed, and therefore where spatial specificity may be tightest. We test this hypothesis with a measurement of the BOLD spatial resolution that uses a spatial activation pattern represented on the flattened cortical surface reconstruction, and determine the fidelity of this pattern as a function of cortical depth. To achieve this measure of the tangential Point-Spread Function (PSF) and its potential changes as a function of cortical depth we designed a novel resolution stimulus to evoke a desired activity pattern within primary visual cortex (V1) that exploits a recent model of the visual topographic mapping, as well as the columnar organization of receptive field centers in V1. This allowed us to use the known visuotopic layout in humans to generate a suitable activation pattern for our measurements. We also implemented a method for generating surface reconstructions at intermediate cortical depths, and utilized a recently developed alignment tool to register the fMRI images to the anatomical images used to create the laminar surfaces. These new analysis tools, combined with 1 mm isotropic voxel sizes and reduced spatial distortion achievable with fast head gradient coils and highly-parallel receive coil arrays—together with the strong signal changes afforded at 7 T field strength—show that targeted sampling of central cortical layers with small gradient-echo EPI voxels can improve the spatial specificity of the fMRI signal, as expected from our knowledge of cortical vascular anatomy. Preliminary accounts of these results have appeared in abstract form [Polimeni et al., 2009a,b].

Methods

Study participants

Six neurologically normal adults (three female) volunteered to participate in the study. Informed consent was obtained from each participant prior to the experiment in accordance with our institution's Human Research Committee. Subjects were compensated financially $50 USD for their participation. Each participant had either normal or contact lens corrected vision.

Visual stimulation

To produce a desired spatial pattern of activation on the V1 surface, a visual stimulus was designed using a model of the visuotopic mapping within V1 to predict the proper stimulus shape in the visual field. We sought a target activity pattern designed to be easily identifiable with recognizable features and with minimal symmetry to reduce ambiguity. Thus the figure “M” was chosen as the target activity pattern. To calculate the visual stimulus that would give rise to this pattern on the cortical surface (see Fig. 1A), we employed a model of the visuotopic mapping given by the known complex-logarithm transformation that describes the topographic mapping of the visual field onto the cortical surface [Schwartz, 1980]. This model is described by two model parameters that control the magnification and local anisotropy of the mapping [Polimeni et al., 2006]. To generate the stimulus, we predicted the mapping in our subjects using this model and assuming the average parameter values (a = 0.7°, γ1 = 0.9) measured in a previous fMRI study across a population of subjects [Polimeni et al., 2005]. We used the inverse of this mapping to generate the stimulus shapes for each visual hemifield that would give rise to the target activity pattern in the contralateral cortical hemisphere. Fig. 1 shows the stimulus conditions. Condition A consisted of the inverse-warped figure filled with a black-and-white spatial noise pattern [Polimeni et al., 2005], superimposed on a gray background. Condition B consisted of the negative image of the same inverse-warped figure filled with background gray superimposed on a flickering spatial noise background (see Fig. 1B). Before each scan, we randomly generated 10 independent spatial noise patterns, and during the stimulus presentation the ensemble of these patterns was presented in sequence, while counterphase flickering at 8 Hz (implemented by presenting each generated pattern for 125 ms, followed immediately by its contrast-reversed version for the same amount of time) to strongly drive V1 neurons. The random spatial noise pattern is advantageous since each pattern it contains energy within small range of spatial scales at each point in the visual field, and the ensemble of patterns contains energy at all spatial orientations at each point in the visual field [Polimeni et al., 2008], thus this stimulus provides a band of spatial scales and all orientations to drive the neurons within each voxel. The luminance of the stimulus ranged from 1.6 to 445 cd/m² such that the luminance contrast modulation of the noise pattern was approximately 99%, and the background gray luminance was approximately 135 cd/m². The stimulus presentation followed a block design schedule where each stimulus condition was presented alternately for 15 s separated by 10 s of background gray. Each functional run consisted of three full stimulus blocks, where one block consisted of the progression: condition A (15 s) — rest (10 s) — condition B (15 s) — rest (10 s). Therefore, BOLD data was acquired for 150 s per run.

Fig. 1

The resolution stimulus. (A) The letter “M” was warped by the inverse of the known complex-logarithm transformation that describes the topographic mapping of the visual field onto the cortical surface in order to evoke the desired activity (more ...)

The stimulus was rear-projected on a screen mounted at the end of the patient table, and the subjects viewed the stimulus through a mirror mounted on the receive coil array. The confined space within the head gradient coil limited the available field of view, and the tight clearance between the receive coil array and the transmit coil in addition to the receive coil interface electronics and housing blocked the bottom circular segment of the lower visual field (less than 1° measured along the vertical meridian). For each subject we could not detect any clipping of the activation pattern in either the single-condition or differential responses, thus no evidence of a reduced field of stimulation was apparent in the data. The stimulated area (projection screen) subtended approximately 24° of visual angle along the horizontal meridian and 14° along the vertical meridian (i.e., extended to 12° horizontal eccentricity and 7° vertical eccentricity). The stimulated area of the “M” pattern subtended approximately 8° of visual angle (i.e., 4° eccentricity). The subtended angles depend slightly on the head size, since this altered the distance from the eye to the mirror, however the variability of the distance of the subject's eye to the projector screen was minor (the maximum difference across subjects was approximately 7%).

A simple fixation task was provided to aid and evaluate subject fixation performance, and to maintain subject vigilance throughout the experiment. The task consisted of a small red dot (approximately 10′ in diameter) located at the center of the visual field that changed luminance randomly over time with a mean interval of 2 s. The subjects were instructed to press a button when the fixation point changed luminance. The subject response was recorded as correct if they responded within 700 ms of the luminance change. During the anatomical scan at the beginning of each session, the subjects practiced the task while performance was continuously monitored. The task difficulty was adjusted by resizing the radius of the fixation dot until the subjects achieved 80% performance. At the end of each scan session the performance score was displayed on the stimulus screen to provide feedback to the subject. In the fMRI analysis, runs with performance below 60% percent were discarded. Time permitting, any discarded runs were repeated at the end of the session after allowing the subjects a short rest within the scanner. Each functional run was limited to 2 min 48 sec to provide the subject with a small rest period to aid in subject fixation performance as well as to limit the lost scan time when a given run was discarded due to either poor fixation performance or excessive head motion.

Anatomical MRI data acquisition

T1-weighted volumetric anatomical data were acquired during each functional (7 T) session for accurate positioning of the fMRI slices using a 1 mm isotropic multi-echo MPRAGE (MEMPRAGE) pulse sequence to limit the susceptibility distortions while maintaining high image SNR [van der Kouwe et al., 2008], and a similar acquisition was collected for each subjected on a separate anatomical session on a 3 T scanner for cortical surface reconstructions. The acquisition parameters for the 3 T data set were as follows: TR=2510 ms; four echoes with TE=1.64 ms, 3.5 ms, 5.36 ms, and 7.22 ms; TI=1200 ms; flip angle = 7°; FOV=256×256×176 mm³, 256×256×176 matrix; bandwidth = 651 Hz/pixel; and R=3 acceleration with online GRAPPA reconstruction using 32 reference lines.

Functional MRI data acquisition

All fMRI data was acquired on a 7 T Siemens whole-body system (Siemens Healthcare, Erlangen, Germany) equipped with AC84 head gradients (80 mT/m maximum gradient strength and 400 T/m/s maximum slew rate). A custom-built 32-channel RF loop coil head array was used [Wiggins et al., 2006] for reception, and a custom-built detunable band-pass birdcage coil for transmit. For the BOLD acquisition we used a standard single-shot gradient-echo interleaved multi-slice EPI protocol: TR/TE/flip=2500 ms/24 ms/90°, fat saturation, FOV=192 mm × 192 mm, a 192 × 192 matrix, bandwidth = 2005 Hz/pixel, trapezoidal readout gradients with 26% ramp sampling, nominal echo spacing of 0.69 ms, and 3-fold acceleration yielding an effective EPI echo spacing (with the GRAPPA acceleration) of 0.23 ms. For each slice we acquired 128 reference lines for parallel imaging reconstruction. The images were reconstructed with the standard online Siemens EPI and GRAPPA reconstruction. Including the initial 10 s block “dummy scans” at the beginning of each scan to allow T1 steady state to be achieved, and the 3 TRs of (three-shot) reference lines acquired immediately following the dummy scans for the GRAPPA reconstruction, the wall-time for each run was exactly 2 min and 47.5 s (or 167.5 s). The 1 mm in-plane resolution images covered the primary visual cortex using 34 1 mm thick oblique-coronal slices positioned orthogonal to the calcarine sulcus with a 20% slice gap. A slice gap was chosen to reduce the potential for slice cross-talk, but results in a larger sample spacing (by 200 μm) in the slice direction; however the voxels themselves are still nominally isotropic. The thin slices were achieved using a modification to the pulse sequence to allow for a longer (7.68 ms filtered sinc) excitation pulse. Conventional B0 field maps derived from phase differences between gradient echo images acquired at TE = 4.22 and 5.24 ms were also acquired at each EPI slice location.

Functional analysis

Prior to statistical analysis, all functional data was motion corrected to the first volume of the first run using MCFLIRT [Jenkinson et al., 2002], then the best-fitting linear trend was removed from each voxel time series. No slice-timing correction was applied. The statistical analysis on the resulting time series data was carried out using FILM with local autocorrelation correction [Woolrich et al., 2001] through modeling both stimulus condition time courses and assuming a canonical hemodynamic response in the form of a Gamma variate with 6 s mean lag for the delay and 3 s standard deviation for the width. The z-statistics for each voxel were computed for the responses to both conditions separately as well as the contrast consisting of the difference between the responses to the two conditions. No spatial smoothing (volume- or surface-based) was applied to the data at any stage of the analysis.

Laminar surface generation

Surface reconstructions of the interface between the cortical gray matter and underlying white matter and the interface between the cortical gray matter and pial surface were automatically generated from the 1-mm MEMPRAGE anatomical data collected at 3 T using FreeSurfer [Dale et al., 1999; Fischl et al., 1999, 2001, 2002; Ségonne et al., 2004, 2005], and cortical thickness maps were derived from these bounding surfaces [Fischl & Dale, 2000]. To estimate the relative positions of the cortical lamina within the gray matter, nine additional surfaces were reconstructed within the cortical gray matter at fixed relative distances between the white and pial surfaces determined from the cortical thickness, i.e., the first intermediate surface was located at 10% of the cortical thickness away from the white matter surface, the second intermediate surface at 20%, and so on. The algorithm for generating the intermediate surfaces followed the procedure of Dale et al. [1999] for computing the pial surface reconstruction from the white matter surface. Specifically, we start with the white matter surface generated by FreeSurfer and then propagate the surface vertices along the line segment connecting the gray/white boundary to the pial surface, while enforcing a curvature smoothness constraint, a self-intersection prevention constraint, and a vertex spacing uniformity constraint [Fischl & Dale, 2000]. The vertex spacing uniformity constraint preserves vertex neighbor distances by penalizing nonuniform stretching of the mesh, and thereby minimizes metric differences between the generated surfaces. Thus for each subject and for each cerebral hemisphere, a collection of 11 surfaces were computed corresponding to the white matter surface, the pial surface, and the nine equally spaced intermediate surfaces. Because each intermediate surface and the pial surface are generated from the initial white matter surface, there is a natural vertex correspondence across all surfaces, and the positions of the corresponding vertices at any location along the white matter surface can be aggregated to form contours through the cortical depths, which can be used to form laminar profiles.

Registration of functional to anatomical surfaces

The laminar surfaces are derived from the MEMPRAGE anatomical data, so it is necessary to align the EPI volumes to these surfaces. Since the surface analysis on these sub-surfaces is the primary goal, we use the corresponding laminar surface information in the EPI images to guide the registration. The principal laminar surface readily identifiable in the EPI volumes is the gray-white surface. We employed a recently proposed boundary-based registration method [Greve & Fischl, 2009] that identifies the interface between gray matter and white matter in the EPI data and then calculates a rigid transformation that registers this interface in the EPI data to the corresponding surface reconstruction from the anatomical data.

Laminar analysis of the functional data

The time-series statistics calculated for each voxel in the functional volume were then transferred onto the collection of surface reconstructions using the transformation produced from the functional-to-anatomical registration by projecting either z statistics or the percent signal change of each voxel intersecting a surface onto the corresponding surface using nearest-neighbor interpolation.

For each subject, a gray matter partial volume map was generated for each hemisphere to quantify the expected contribution of gray matter signal relative to other white matter or CSF/pial contributions. The map was generated from the pial and white-matter surfaces by filling the space between and around these surfaces with 0.25×0.25×0.25 mm voxels labeled as white, gray or CSF. This high resolution model was then resampled with trilinear interpolation into a volume with 1 mm voxels identical to the BOLD acquisition aligned to the anatomical space with the computed registration. This produced three volumes with the same geometry as the fMRI data containing the relative contribution of each compartment such that at any given voxel the sum of the corresponding voxels from the three volumes was 100%. Then the partial volume map for each component (gray, white, and CSF) was projected onto each of the cortical surface reconstructions in the same manner as the functional data. This enabled quantification of the relative contribution of each tissue compartment as a function of laminar depth for each surface vertex. ROI analysis also provided statistics of the mean and standard deviation of the relative contributions limited to the region considered in the functional analysis.

ROI analysis

To quantify the distribution of percent signal change and the temporal features of the BOLD response with V1 as a function of cortical depth, we estimated the position of the full extent of primary visual cortex and its boundaries along the cortical surface using a new atlas method that predicts the location of V1 based on the cortical folding pattern [Hinds et al., 2009]. Because we were able to stimulate only 20° of the visual field, approximately half of the visual cortex was activated by the stimulus, thus we chose a region of interest (ROI) consisting of the predicted V1 beginning at the occipital pole and extending halfway through the predicted region along the axis of the calcarine sulcus in the anterior direction. Aside from the results that pool data within this V1 ROI, no smoothing (volume-based or surface-based) was applied to the data at any stage of the analysis. Furthermore, no spatial normalization was employed, and all analyses were conducted in individual subjects.

For each subject, a patch of the cortical surface consisting of the V1 ROI intersected with the acquisition volume was identified then extracted from each surface in the collection and computationally flattened using a near-isometric flattening algorithm [Balasubramanian et al., 2010] that minimizes flattening distortion. This algorithm provides a per-vertex error measure to assess the flattening accuracy for each surface. In summary, the near-isometric flattening algorithm begins by computing the shortest-path distances between each pair of vertices along the cortical surface mesh (along the faces of the mesh, i.e., the paths are not constrained to lie along the edges of the triangulation). Note that these distances correspond to minimal-geodesic distances when the surface is geodesically convex. The flattening algorithm then computes an optimal configuration of the vertices in the plane such that the two-dimensional Euclidean distances in the plane best match the shortest-path distances measured along the folded surface. Since these surfaces have non-zero Gaussian curvature, it is not possible to flatten into the plane isometrically, but we can characterize the flattening “distortion” by comparing how well the Euclidean distances in the plane of the final flattened surface match the shortest-path distances along the original surface. The metric distortion we report is therefore the relative discrepancy between the shortest-path distances and the Euclidean distances (i.e., |D_original–D_flattened|/D_original), averaged over the flattened patch.

To quantify the spatial fidelity of the measured activation pattern relative to the desired activation pattern as a function of cortical depth, the normalized correlation between the ideal activation pattern and the measured pattern (resulting from the differential response to the two stimulus conditions) was computed in the flattened representation of the cortical surface. Prior to the correlation analysis, the activation template was aligned to the differential activation pattern within V1. So as to not bias the alignment to any particular depth, for each hemisphere we aligned the template to a single activation pattern consisting of the sum of the z-statistics collapsed across all depths, and used this one alignment to compute the spatial correlation at each depth. The alignment was calculated using a multi-scale search algorithm that found the best translation, rotation, and scaling parameters using a phase-only template matching algorithm [Horner & Gianino, 1984]. After the alignment was computed for each hemisphere, we visually inspected the resulting alignment to verify that the alignment process converged and that it produced an accurate alignment of the template to the data in order to prevent a failed convergence from biasing the results of the subsequent correlation analysis.

Results

EPI acquisition

Example 1 mm single-shot gradient-echo EPI images are shown in Fig. 2A. For this subject, the mean off-resonance frequency shift within the V1 gray matter as determined from the field maps was 23 ± 23 Hz across our population of subjects, and the full-width-at-half-maximum (FWHM) of the best-fitting Gaussian to the off-resonance distribution was 36 ± 19 Hz. For the 0.23 ms effective echo spacing this corresponded to an average pixel displacement of 0.36 mm ± 0.18 mm across all subjects.

Fig. 2

(A) Example EPI slices from a BOLD fMRI scan. (B) Registration of 1-mm isotropic MEMPRAGE anatomical data with EPI data. White matter surface reconstruction is superimposed as yellow contour.

Registration and laminar surface generation

To demonstrate the accuracy of the boundary-based registration algorithm for aligning the EPI slices into the coordinates of the anatomical data, the surface reconstructions of the interface between the gray and white matter computed from the MEMPRAGE data are overlaid on top of the individual BOLD-weighted images in Fig. 2B. The surface reconstruction intersects the gray-white interface detected in the EPI slices.

A family of 11 surface representations was generated for each cerebral hemisphere. Fig. 3A shows a subset of these surfaces confined to the cortical gray matter superimposed on a slice of the MEMPRAGE anatomical data. The surfaces are evenly spaced throughout the cortical thickness and extend throughout V1. Figures 3B and C show an example of how the cortical surface models are positioned in a relatively flat section of cortex within the visual cortex. The method evenly spaces the surfaces across the cortical thickness, and therefore provides a uniform sampling of cortical depth. Due to the folded nature of the cortex, each intermediate surface contained differing amounts of curvature and total surface area throughout the V1 ROI. The lowest absolute value of the mean curvature within V1 was found for the surface at depth 6—corresponding to the approximate location of Layer IV—as has been previously observed [Van Essen & Maunsell, 1980]. The minimum of the absolute value of the curvature is expected in the center-most surfaces since this location is maximally distant from the high curvature pial and white matter boundaries.

Fig. 3

Laminar analysis of the fMRI activation. (A) Given the cortical surface models, the activation patterns can be projected onto intermediate surfaces lying within the thickness of the cortical gray matter. With small fMRI voxel size, many voxels lie in (more ...)

Because the average cortical thickness within V1 for this population of subjects was 1.8 mm ± 0.1 mm and the EPI voxel size was 1 mm, the nine intermediate surfaces sampled the cortical depths more finely than the functional acquisition. However, due to the folded nature of the cortex, the number of EPI voxels intersecting the laminar surfaces varied with depth. Because the curvature in V1 is predominantly negative there were consistently fewer voxels (on average 1% across subjects within V1) intersecting the white matter surface than the pial surface, and this albeit small difference was significant in our population (p<0.05, one-tailed).

Spatial fidelity of resolution pattern across cortical depths

The statistical analysis of the differential responses to the two conditions of the resolution stimulus is shown in Fig. 4 for one subject as a z-statistic map sampled on intermediate cortical surface reconstructions across five cortical depths. The activity becomes stronger for pixels nearer to the pial surface, with the maximum z statistics near the pial surface and the minimum near the white matter surface. For this subject, the average z statistic value within the V1 ROI dropped by 45% from the pial surface to the white matter surface (36±9% across the population). The fidelity of the resolution pattern, however, degrades near the pial surface. Together, these findings are consistent with reduced activity from partial volume interactions with white matter, and the dominance of macrovasculature at the pial surface. The pattern is maximally coherent in the deeper laminae where the overall activity strength is also minimal. This suggests that the spatial specificity is lowest near the pial surface where the extra-vascular signal is dominated by the responses from within large draining veins.

Fig. 4

The activity pattern of the “M” stimulus calculated from two 2.5 min. acquisitions (i.e., 5 min. total acquisition time) resulting from differential imaging of the two stimulus conditions, and displayed on the inflated cortical surface. (more ...)

Fig. 5A shows an example flattened patch overlaid with vertex position at which the z statistic exceeded a threshold of 2.0. The corresponding map of the flattening error is shown in Fig. 5B. The metric distortion assessed on the flattened patch using the per-vertex flattening error measure calculated from the near-isometric flattening [Balasubramanian et al., 2010] and was found to be on average 3.5% (± 0.7% std.) and everywhere less than 6.8% within the V1 ROI. After computing the optimal alignment between the positions of the activity pattern template and the measured differential activation pattern on the flattened cortical surface as described above, the correlation between the template and the differential activation pattern was found to decrease consistently with depth. The normalized correlation averaged over our population of six subjects (i.e., 12 hemispheres) is shown in Fig. 5C. The correlation between the template and the activation pattern decreased consistently with depth, with a 35% increase from the pial surface to the white matter. Thus the fidelity of the pattern decreased consistently as the pattern was sampled closer to the pial surface, confirming the trend that is visually apparent in Fig. 4. Fig. 5D shows the standard deviation of the normalized correlation across our population of subjects, demonstrating that the deterioration of the spatial fidelity with proximity to the pial surface varied across subjects and that this variability was more pronounced for sampling depths close to the pial surface.

Fig. 5

Estimation of spatial spread of activity across cortical depths. (A) Example near-isometric flattening of a patch of the cortical surface with the above-threshold activated vertices labeled, and (B) the corresponding flattening error map. The metric distortion (more ...)

To determine the effect of cortical depth on the magnitude of the BOLD signal changes, we computed the average percent signal change and the spatial variability of the percent signal change in the responses to each stimulus condition over the V1 ROI as a function of cortical depth. Fig. 6A demonstrates that the average percent signal change within the ROI increases monotonically from the deep layers up to the shallow layers and is maximal at the pial surface, consistent with a peak response strength stemming from the large draining vessels on the cortical surface. Therefore a trade-off exists across depths between the spatial accuracy and the response strength. In Fig. 6B, the spatial variability of the percent signal change (i.e., the degree to which the percent signal change varies with position along the cortex) is plotted across cortical depths and the trend is the same. This provides an indirect measure of the degradation of spatial fidelity of the activity pattern (given that the response strength is more spatially inhomogeneous—and therefore more erratic—near the pial surface). Thus both the strength and spatial variability of the signal changes are maximal at the pial surface, consistent with a situation in which large signal changes are present near surface vessels and moderate changes are present between vessels, which may contribute to the high spatial variability in BOLD signal strength at the pial surface.

Fig. 6

(A) The percent signal change across depths averaged over all vertices within the V1 ROI. Peak percent signal change appears near the pial surface, in agreement with the larger z-statistics seen in Fig. 4 near the surface vessels, despite the lower spatial (more ...)

Partial volume effects as a function of cortical depth

For each subject, a partial volume map was generated to quantify the relative contribution of gray and white matter tissues and CSF in each laminar analysis group. The average contributions of gray matter, white matter, and CSF to the voxels included in the functional analysis over all subjects are presented in Fig. 7 as functions of cortical depth. This analysis demonstrates that, for our 1 mm isotropic voxel acquisition, only about 35% of the EPI voxels are contained entirely within the cortical gray matter. This indicates that both the breakdown of the spatial fidelity of the activity pattern and the steady decrease in percent signal change with cortical depth may be in part attributable to partial volume effects and the contributions of signals from outside of the gray matter tissue.

Fig. 7

Partial volume effects across depths within a V1 ROI. Average percentage of voxel volume contributed from the three categories: white matter (WM), gray matter (GM), and cerebrospinal fluid (CSF). Only central cortical layers contain gray matter when using (more ...)

To explore the relationship between functional activation significance and the tissue mixture in the voxel, we computed the average functional CNR (in the form of the z statistic) as a function of the relative contribution of each tissue type to the EPI voxels. The results are given in Fig. 8, which indicate that the functional CNR increases with increasing contributions of gray matter to the voxel composition, but only the largest z-statistic values were found in voxels containing a combination of gray matter and CSF, consistent with strong pial vessel contributions for these voxels. However, in both the individual subject and across the population, only slightly smaller z-statistic values are found in voxels that are entirely within the gray matter—voxels where the measured z-statistic was 80% of the maximum contained more than 90% gray matter. Limiting analysis to voxels exclusively containing gray matter therefore results in only a small penalty in effect size, but can improve spatial specificity.

Fig. 8

Functional CNR contributions from gray matter, white matter, and CSF. (A) The relative contributions of the three tissue classes, cortical gray matter, white matter, and surrounding CSF, are plotted against the average z-statistic—representing (more ...)

Single-condition vs. differential imaging maps across layers

The activation pattern presented in Fig. 4 is implicitly generated from the contrast between two stimulus conditions that ideally activate spatially distinct regions of cortex, and is therefore a form of differential imaging. Differential imaging is often employed to enhance the localization of fine-scale spatial features relative to what is achievable with standard single-condition imaging. The two single-condition activation patterns are presented in Fig. 9A in the same format as in Fig. 4 and correspond to the same subject. While the desired pattern of the resolution stimulus can be detected in the single-condition maps, several artifacts common between these two maps are removed when the two single-condition maps are subtracted.

Fig. 9

Activation maps from single-condition responses. (A) The responses to “stimulus condition A” (top) and “stimulus condition B” (bottom), across several cortical depths. (B) The normalized correlation between the desired (more ...)

To compare the relative degradation of the single-condition responses across cortical depths to that of the differential responses, we performed the same spatial correlation analysis between a template of the ideal activity pattern and the measured single-condition responses to the corresponding stimulus. To enable direct comparison with the correlation analysis performed with the differential-condition activation patterns, we used the same computed alignment for the single-condition correlation analysis. The results are provided in Figs. 9B and C across our population of subjects. The normalized correlation was lower than for the differential response, but relatively constant across all cortical depths compared with the analysis of the differential condition maps (shown in Fig. 5).

Discussion

Visuotopic mapping in V1

In this study, knowledge of the visuotopic mapping within cortical area V1 enabled us to impose a desired spatial pattern of activation along the cortical surface without requiring any visuotopic mapping measurements in this group of participants. While the V1 model adopted here assumed model parameter values taken as the average values across a population of subjects, naturally some moderate variability exists [Polimeni et al., 2005]. Thus any mismatch between the model and the true mapping within an individual can corrupt the fidelity of the desired activity pattern and lead to its distortion. Furthermore, the inverse-warping of the stimulus assumed that the subjects were fixating the exact center of the visual field, therefore any errors in eye position such as fixation offset or eye movements are a confound for the experiment. Typical patterns of eye movements have been shown to be idiosyncratic across subjects—both fixational eye movements [Ditchburn, 1973] as well as the systematic eye fixation errors that occur even when the subject perceives that they are staring at the target and the eye is stationary [Putnam et al., 2005]. In the case of either a difference in a subject's retinotopic mapping from the population average or a systematic and static error in subject eye fixation, the test pattern activation would appear warped when viewed on the flattened cortex.

In addition to serving as a tool for aiding in the characterization of the biological point-spread function of fMRI and the spatial specificity afforded by our acquisition and analysis methods, the success of this technique partially validates the parametric model of V1 topography used to design the stimulus. The accurate prediction of the mapping in novel subjects was enabled by the consistent structure of the topographic mapping across individuals [Polimeni et al., 2005]. Recent demonstrations of the regularity of the shape of V1 [Hinds et al., 2008] and of the consistent position of V1 relative to the cortical folding pattern [Fischl et al., 2008; Hinds et al., 2009] show that a remarkable regularity of V1 structure and function exists in the geometry of this primary sensory area.

Differential imaging for feature localization

Often in neuroimaging, differential imaging is employed to localize features in topographic maps [Bonhoeffer & Grinvald, 1993; Yacoub et al., 2007]. These techniques operate by presenting pairs of stimuli that are known to preferentially activate spatially distinct populations of neurons, and thus the corresponding response patterns are subtracted in order to eliminate any common responses that are present under both stimulus conditions. The success of this method depends on prior knowledge of the existence of a topographic mapping, and also requires that the responses of the underlying neural activation possesses spatial linearity [Hansen et al., 2006]. If, for example, the point-spread function of the response is non-zero-mean, as is the case for activations measured near the cortical surface near draining veins, it is possible that the assumptions are invalid, and therefore differential imaging can introduce spatial bias. Regardless, subtracting different conditions can aid in the spatial localization of features of the topographic mapping by leveraging the one-to-one correspondence between visual field and position in V1, but the improvements in spatial fidelity come from the validity of this prior information.

The analysis of differential responses to visual stimuli presented in this study was in part motivated by the common use of similar differential contrast in fMRI in the visual system, and in particular in high-resolution studies [e.g., Yacoub et al., 2007]. Here we tested how well differential imaging can remove the common responses between two stimulus conditions to improve spatial localization, and compared the performance of differential imaging across several cortical depths. The comparison between single-condition and differential activation maps presented in Fig. 9 demonstrates that indeed some of the spatial spread in the response is removed in the differential map. However, the spatial fidelity of the differential response sampled near to the pial surface vessels, as seen in Fig. 4, exhibited increased spatial errors compared to the responses sampled far from the pial surface vessels. Thus, differential imaging cannot fully remove the large, systematic errors that corrupt data sampled near to the pial surface.

Surface-based laminar analysis

In this study, the combination of small voxel sizes in the fMRI acquisition and accurate boundary-based registration enabled sampling the BOLD signal at multiple cortical depths. By generating intermediate surface reconstructions between the boundaries of the cortical gray matter, a tangential sectioning of the functional data could be visualized and analyzed independently for multiple depths. The positions of these intermediate surfaces were chosen based on the normalized distance between the white matter boundary and pial surface boundary (subject to mild smoothness constraints). Unfortunately, these intermediate surfaces are unlikely to align perfectly with the cortical lamina, since the laminar features are known to vary in depth as a function of the local curvature of the cortex. Thus, while the cortical layers are known to lie at fixed depths within naturally flat sections of cortex within a given cortical area, the relative thickness of the supragranular layers and infragranular layers varies over the cortex and is strongly correlated with the features of the folding pattern—that upper layers are consistently thinner within the crests of cortical gyri, and the lower layers are thinner within the fundi of cortical sulci has been demonstrated in both macaque [Van Essen & Maunsell, 1980; Hilgetag & Barbas, 2006] and human [Fatterpekar et al., 2002]. Cortical thickness also varies consistently across folding pattern: it is thickest on crowns of gyri, thinnest in fundi of sulci, with an intermediate thickness along the banks of sulci [Fischl & Dale, 2000; Welker, 1990]. This suggests that the next step in our laminar analysis stream should be to apply a correction factor, which is a function of the local cortical curvature, to the placement of the laminar surfaces so that the surfaces better approximate the position of the histological laminar features rather than simply being determined from fixed percentages of the cortical depth.

Another method that has been proposed for predicting the cortical laminar morphology from the boundary surfaces calculates solutions to the Laplace equation—or, alternatively, repeatedly applies the corresponding Green's function, the heat kernel—by placing boundary conditions on the gray matter and pial surface borders, and assigns the level lines or equipotentials of the solution to the cortical layers [Jones et al., 2000]. This method has the additional advantage of supplying streamlines of the solution that are a family of contours orthogonal to the potential lines, which do not intersect and thus are convenient to use to generate laminar profiles. While intuitively the equipotential contours are expected to bunch where the boundaries exhibit high curvature, causing a local laminar compression, in practice this approach does not accurately capture the laminar morphology [Annese et al., 2004]. More work is needed to infer the positions of the cortical lamina within the gray matter. Once accurate predictions of laminar morphology can be made from the curvature of the cortical pattern, it will be possible to simulate tangential sectioning and tissue flattening that has been employed in classical anatomical studies that revealed large-scale organization of the visual cortex, which could not be fully appreciated either through the point-sampling afforded by micro-electrode physiological recordings or though histology on serial tissue sections [e.g., Levay et al., 1975; Tootell et al., 1982; Horton & Hedley-Whyte, 1984].

In addition to providing an understanding of the sources of resolution degradation in fMRI, the readily recognizable retinotopic pattern also validated the overall ability of the fMRI data to accurately align to a cortical surface model, and that the many steps involved in a high-resolution laminar analysis across the folded cortical surface can be accurately performed. For example, if the local alignment between the functional data and the surface models systematically varied across the test pattern by a substantial fraction of the voxel size, the regions of the pattern that were misaligned would show substantially reduced activation for the surfaces near the white matter and pial boundary where misalignments would cause non-activated WM or CSF voxels to be identified as gray matter voxels at these boundary lamina.

Vascular density closely follows the laminar morphology in both macaques [Zheng et al., 1990; Weber et al., 2008] and humans [Duvernoy, 1981; Lauwers et al., 2008], with the band of highest density found at cortical Layer IV. In area V1 in particular, it has been claimed that Layer IV is widest [Blinkov & Glezer, 1968], and is the site of the highest vascular density [Duvernoy, 1981], of all cortical areas. As with the thickness of the upper and deep cortical lamina, significant variations in vascular density are correlated with cortical folding pattern, with higher vascular density appearing at the crowns of gyri and lower vascular density at the fundi of sulci [Duvernoy, 1981; Lauwers et al., 2008]. Thus the spatial specificity of the BOLD signal is likely to track the variations in vascular density and will depend on the position within the folding pattern as well as the depth within the cortical gray matter.

Our current inability to sample precisely from Layer IV suggests that we are sampling from cortical compartments with varying vascular density, thus the spatial accuracy of our laminar assignments are likely to vary slightly with position. Also, the 1 mm isotropic voxels employed in this study may not be sufficiently small to allow for measuring these small features of the BOLD response across layers in the primary visual cortex, which is thinner compared to many cortical areas [Blinkov & Glezer, 1968]. Taken together, these issues may explain why we do not see evidence of a local peak in the BOLD response at Layer IV, as seen in other studies where locally flat areas were analyzed [Ress et al., 2007; Koopmans et al., 2009]. We are currently investigating extensions of our laminar-surface-based methods to directly sampling from specific cortical layers by modeling the relationship between local curvature and laminar depth. Together with smaller voxel sizes, this could potentially enable improved detection of specific vascular layers.

Small-voxel fMRI acquisition at 7 T

While gradient-echo EPI is the most commonly employed pulse sequence for fMRI, spin-echo-based sequences have been shown to be preferentially sensitive to extra-vascular signal changes near smaller vessels and capillaries [Boxerman et al., 1995], and are therefore advocated for high-resolution fMRI at high field strengths [e.g., Duong et al., 2003; Yacoub et al., 2007]. Unfortunately there are several disadvantages to spin-echo EPI, including its reduced overall sensitivity compared to gradient-echo EPI, and the required increase in power deposition and Specific Absorption Rate which must be limited for subject safety. Furthermore, dielectric wavelength shortening effects at high field strengths present challenges for the spin-echo sequence resulting in a loss of signal for spatial positions where the B₁ transmit inhomogeneity causes reductions from the prescribed flip angles. Additionally, spin-echo EPI does not provide ideal spin-echo-based T₂ contrast for all spatial frequency components within the image since the EPI readout requires a substantial period of time compared to T₂*, thus causing T₂* decay at the edges of k-space. Thus, while the center of k-space is fully refocused and is purely T₂ weighted, the edges of k-space, containing important information on fine (high-spatial-frequency) details in the image possess a combination of T₂ and T₂* weighting [Duyn, 2004; Keilholz et al., 2005; Goense & Logothetis, 2006; Poser & Norris, 2007]. Our use of small gradient-echo EPI voxels and anatomically-informed sampling to restrict the analysis to BOLD responses within the central and deep cortical layers away from the large surface vessels provides a means to selectively sample BOLD responses originating from a smaller distribution of vessel diameters compared to the distribution encountered when the pial surface contaminates the cortical voxels through partial volume effects. Although this approach is not able to fully exclude all non-capillary sized vessels (since the frequent diving venules cannot be excluded), the method yields improved localization of functional boundaries and thus provides some of the spatial specificity achieved by the spin-echo EPI method without the associated losses of sensitivity.

Our work shows that the effective biological spatial resolution of fMRI can be improved by sampling within the both the deep and the central layers of cortex, and thereby explicitly excluding the pial vessels, and that the central layers near to the capillary bed in Layer IV exhibit high spatial specificity while retaining high CNR. This is consistent with previous animal studies that reported both a degraded gradient-echo BOLD spatial fidelity in the vicinity of the pial surface and a three-fold stronger gradient-echo response compared to spin-echo BOLD, and suggested that gradient-echo BOLD may be the preferred method for mapping small spatial features in human imaging in parenchymal regions excluding the large pial vessels [Moon et al., 2007]. While it may be possible that the spatial specificity of BOLD in Layer IV is tighter than in other layers due to the higher capillary density (and thus the smaller spacing between vessels), our results do not demonstrate improved specificity in Layer IV compared to deeper layers—the specificity near Layer IV is only demonstrated to be tighter than more superficial layers. The small improvement in spatial specificity in the deepest layers near the white matter surface may derive solely from the increased distance from the pial vessels. Because the magnetic field disturbance surrounding a long, straight cylinder falls off inversely with squared distance [Ogawa et al., 1993; Boxerman et al., 1995], the largest pial vessels are expected to have an influence over all layers and for the influence to diminish with distance.

Achieving such a fine in-plane voxel size over an image field-of-view the size of a human head requires many time-consuming image encoding steps. Conventional, unaccelerated spatial encoded EPI is thus impractical at 7 T due to the rapid signal decay and increased vulnerability to susceptibility-induced geometric distortions. To ameliorate these factors our acquisitions made use of a custom-built 32-channel receive-only head coil array, which enables acceleration with parallel imaging techniques to reduce the EPI readout times and avoid the severe geometric distortions common in high-field fMRI [e.g., de Zwart et al., 2002]. In addition to reducing susceptibility-induced distortion through R=3 accelerated imaging, we also utilized fast head gradients capable of an EPI echo spacing almost half that of conventional whole-body gradients. The thin slices needed for small voxel volumes helps to reduce through-plane spin dephasing [e.g., Merboldt et al., 2000], which can lead to signal drop-outs in regions in the head near large susceptibility gradients. Although multi-shot techniques, such as segmented EPI, are available for achieving small voxel sizes, they lead to increased scan time due to their lower temporal sampling rate and are more vulnerable to instabilities such as subject motion and physiological noise [e.g., Glover & Lee, 1995].

In the present study, a 20% slice gap was used in the acquisition to reduce spin-history effects that can be more severe at high field strength due to the prolonged gray matter T₁—but at the expense of inducing a coarser spatial sampling (by 200 μm) in the slice direction than the in-plane sampling. Although the voxels are nominally isotropic, and this asymmetry in the spatial sampling may be negligible in the present study, a contiguous acquisition may still be advantageous for accurate localization of activation. In future work we will quantify this trade-off between slice interference and contiguous sampling to determine the optimal sampling strategy.

Improving and characterizing spatial accuracy of fMRI

Although most vessels in the central cortical layers are capillaries oriented roughly randomly with respect to the cortical laminae, diving venules, oriented orthogonal to the cortical laminae, are irregularly distributed and drain deoxygenated blood from the capillaries back to the surface [Duvernoy et al., 1981]. These vessels penetrate Layer IV and pool deoxygenated blood from a large region, and can therefore potentially reduce radial spatial resolution in gradient-echo EPI. The BOLD contamination from the largest of these vessels can be suppressed by identifying voxels containing diving venules as dark spots in the images and removing those voxels from subsequent analysis [Yacoub et al., 2001a; Barth & Norris, 2007]. Thus small voxels centered in Layer IV with vein suppression will at least be partially shielded from partial volume effects of larger vessels.

While techniques such as the “initial dip”, flow-based functional contrast, and spin-echo EPI have been suggested for improvement of the biological point-spread function of BOLD, our study suggests that gradient-echo EPI with small voxels limited to the center lamina of the cortex gray matter can also significantly improve spatial localization. Furthermore, the “vein avoidance” strategy of limiting the study to central lamina can be combined with these other methods for resolution enhancement. Future experiments will seek to systematically identify and correct errors apparent in the spatial pattern, then correct these errors and observe the changes in the fidelity of the spatial pattern of activity.

Typically the BOLD PSF refers to the ability of the fMRI signal to localize activity tangential to the cortical surface, which we will refer to as the tangential point-spread-function PSF_tang, whereas with the availability of laminar analysis the characterization of the point-spread function radial to the cortical lamination, PSF_radial, becomes possible. These two measures of point-spread are possibly distinct. For example, radial veins within the cortex would not be expected to significantly contribute to PSF_tang, instead corrupting only PSF_radial. And, as this study shows, the presence of tangential vessels of relatively large size at the pial surface leads to a corruption of PSF_tang which is a function of the laminar depth. Thus, the BOLD point-spread in a cortical area is likely not well parameterized by even two PSFs (PSF_tang and PSF_radial), but by these two measures as a function of cortical depth.

Previous studies [e.g., Engel et al., 1997; Menon & Goodyear, 1999; Leite et al., 2004; Shmuel et al., 2007] have reported measures of PSF_tang taken by inducing spatial activation patterns in the visual cortex. The most recent measurement of spatial resolution made by Shmuel et al. [2007] using gradient echo imaging at 7 T resulted in a reported PSF FWHM of 2.34 mm (±0.20 mm). However, given the large systematic displacements of activity that occur when sampling near to the large pial vessels [Olman et al., 2007], modeling the BOLD point-spread with a single zero-mean kernel, and summarizing resolution with its FWHM, may not provide a complete characterization of spatial specificity of the BOLD response given these non-zero-mean errors. In this study, our use of a resolution pattern imposed on the cortical surface enabled an analysis of the spatial accuracy of the response through correlating the measured pattern with a template of the ideal response. With this correlation analysis, we were able to assess the spatial accuracy of the BOLD response as a function of cortical depth and to demonstrate the relative degradation in spatial fidelity of this pattern as it was sampled closer to the pial surface. These results provide both a clear demonstration of how veins disrupt the assignment of tangential functional boundaries, and a measure of how this process is improved by vein avoidance.

Additionally, given the highly asymmetric, stereotyped organization of the vascular density tangential and radial to the cortical surface, it is natural to ask whether the local vascular regulation and blood delivery density is spatially specific enough to allow for fMRI resolution of laminar level neuronal signals. While this study provides a methodological foundation for future exploration of this topic, it does not directly address this issue.

Conclusion

Here we have demonstrated the use of a novel laminar surface-based analysis method over the folded cortex and apply it to measuring the spatial accuracy of high-resolution fMRI using a gradient-echo EPI acquisition at 7 T. While techniques such as flow-based functional contrast, spin-echo EPI at high field strengths, and the “initial dip” or early components of the measured hemodynamic response have all been suggested for improvement of the biological point-spread function for fMRI measurements, our study suggests that gradient-echo EPI with small voxels limited to the central and deep lamina of the cortical gray matter can also significantly improve spatial localization.

Previous studies of laminar differences in the fMRI signal have used hand-drawn profiles through cortical depths, orthogonal to the gray matter boundaries in small, flat patches of cortex away from the cortical folds [Ress et al., 2007; Koopmans et al., 2009]. Our novel surface-based laminar analysis method enabled sampling the fMRI signals within several intermediate cortical layers to quantify spatial specificity as a function of cortical depth or proximity to the draining veins. This surface-based laminar analysis is capable of analyzing patterns of activity within the cortex at a particular cortical depth over the full extent of the cortex, and is therefore not limited to a small region of interest or restricted to a flat patch of cortex exhibiting little or no curvature.

The new surface-based analysis techniques enable the accurate localization of small voxels within the central layers of the cortical gray matter, which we have shown is required to achieve high spatial specificity. For this demonstration, we exploited the known spatial structure of the visuotopic mapping in V1 to impose a novel resolution stimulus, which provided a clear characterization of the spatial accuracy across cortical depths. Due to our knowledge of its topographic structure, the visuotopic mapping in V1 therefore can provide a valuable testing ground for measurements of spatial resolution in fMRI.

Acknowledgments

We would like to thank Dr. Thomas Benner for his assistance with the EPI protocol, Ms. Kyoko Fujimoto for her help conducting the experiments, Dr. Graham Wiggins for the 32-channel array coil, and Mr. Thomas Witzel and Drs. Oliver Hinds, Mukund Balasubramanian, and Eric Schwartz for useful discussions. This work was supported by the Athinoula A. Martinos Center for Biomedical Imaging, the NIH NCRR grant P41 RR14075. Further support was provided by the National Institute for Biomedical Imaging and Bioengineering (R01 EB006758, and R01 EB006847), the National Institute on Aging (AG02238), the National Institute for Neurological Disorders and Stroke (R01 NS052585-01) and by the Autism & Dyslexia Project funded by the Ellison Medical Foundation.

Footnotes

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

References

Annese J, Pitiot A, Dinov ID, Toga AW. A myeloarchitectonic method for the structural classification of cortical areas. NeuroImage. 2004;21(1):15–26. [PubMed]
Balasubramanian M, Polimeni JR, Schwartz EL. Near-isometric flattening of brain surfaces. NeuroImage. 2010 in press. [PMC free article] [PubMed]
Bandettini PA, Wong EC, Jesmanowicz A, Hinks RS, Hyde JS. Spin-echo and gradient-echo EPI of human brain activation using BOLD contrast: a comparative study at 1.5 T. NMR Biomed. 1994;7(1–2):12–20. [PubMed]
Barth M, Norris DG. Very high-resolution three-dimensional functional MRI of the human visual cortex with elimination of large venous vessels. NMR Biomed. 2007;20(5):477–84. [PubMed]
Blinkov SM, Glezer II. The Human Brain in Figures and Tables. Basic Books, Inc., Plenum Press; New York: 1968.
Bonhoeffer T, Grinvald A. The layout of iso-orientation domains in area 18 of cat visual cortex: optical imaging reveals a pinwheel-like organization. J Neurosci. 1993;13(10):4157–80. [PubMed]
Boxerman JL, Hamberg LM, Rosen BR, Weisskoff RM. MR contrast due to intravascular magnetic susceptibility perturbations. Magn Reson Med. 1995;34(4):555–66. [PubMed]
Dale AM, Fischl B, Sereno MI. Cortical surface-based analysis. I. Segmentation and surface reconstruction. NeuroImage. 1999;9(2):179–94. [PubMed]
de Zwart JA, van Gelderen P, Kellman P, Duyn JH. Application of sensitivity-encoded echo-planar imaging for blood oxygen level-dependent functional brain imaging. Magn Reson Med. 2002;48(6):1011–20. [PubMed]
Ditchburn RW. Eye-movements and Visual Perception. Oxford University Press; 1973.
Duong TQ, Yacoub E, Adriany G, Hu X, Ugurbil K, Kim SG. Microvascular BOLD contribution at 4 and 7 T in the human brain: gradient-echo and spin-echo fMRI with suppression of blood effects. Magn Reson Med. 2003;49(6):1019–27. [PubMed]
Duvernoy HM, Delon S, Vannson JL. Cortical blood vessels of the human brain. Brain Res Bull. 1981;7(5):519–79. [PubMed]
Duyn J. Specificity of high-resolution BOLD and CBF fMRI at 7 T. Magn Reson Med. 2004;51(3):644–5. [PubMed]
Engel SA, Glover GH, Wandell BA. Retinotopic organization in human visual cortex and the spatial precision of functional MRI. Cereb Cortex. 1997;7(2):181–92. [PubMed]
Ernst T, Hennig J. Observation of a fast response in functional MR. Magn Reson Med. 1994;32(1):146–9. [PubMed]
Fatterpekar GM, Delman BN, Boonn WW, Gultekin SH, Fayad ZA, Hoff PR, Naidich TP. MR microscopy of normal human brain. Magn Reson Imaging Clin N Am. 2003;11(4):641–53. [PubMed]
Fischl B, Dale AM. Measuring the thickness of the human cerebral cortex from magnetic resonance images. Proc Natl Acad Sci U S A. 2000;97(20):11050–5. [PubMed]
Fischl B, Liu A, Dale AM. Automated manifold surgery: constructing geometrically accurate and topologically correct models of the human cerebral cortex. IEEE Trans Med Imaging. 2001;20(1):70–80. [PubMed]
Fischl B, Rajendran N, Busa E, Augustinack J, Hinds O, Yeo BT, Mohlberg H, Amunts K, Zilles K. Cortical folding patterns and predicting cytoarchitecture. Cereb Cortex. 2008;18(8):1973–80. [PubMed]
Fischl B, Salat DH, Busa E, Albert M, Dieterich M, Haselgrove C, van der Kouwe A, Killiany R, Kennedy D, Klaveness S, Montillo A, Makris N, Rosen B, Dale AM. Whole brain segmentation: automated labeling of neuroanatomical structures in the human brain. Neuron. 2002;33(3):341–55. [PubMed]
Fischl B, Sereno MI, Dale AM. Cortical surface-based analysis. II: Inflation, flattening, and a surface-based coordinate system. NeuroImage. 1999;9(2):195–207. [PubMed]
Gati JS, Menon RS. High resolution fMRI using 3D-EPI [Abstract] Proc ISMRM. 2002;10:124.
Glover GH, Lee AT. Motion artifacts in fMRI: comparison of 2DFT with PR and spiral scan methods. Magn Reson Med. 1995;33(5):624–35. [PubMed]
Goense JB, Logothetis NK. Laminar specificity in monkey V1 using high-resolution SE-fMRI. Magn Reson Imaging. 2006;24(4):381–92. [PubMed]
Greve DN, Fischl B. Accurate and robust brain image alignment using boundary-based registration. NeuroImage. 2009;48(1):63–72. [PMC free article] [PubMed]
Hansen KA, David SV, Gallant JL. Parametric reverse correlation reveals spatial linearity of retinotopic human V1 BOLD response. NeuroImage. 2004;23(1):233–41. [PubMed]
Harel N, Lin J, Moeller S, Ugurbil K, Yacoub E. Combined imaging-histological study of cortical laminar specificity of fMRI signals. NeuroImage. 2006;29(3):879–87. [PubMed]
Hilgetag CC, Barbas H. Role of mechanical factors in the morphology of the primate cerebral cortex. PLoS Comput Biol. 2006;2(3):e22. [PubMed]
Hinds O, Polimeni JR, Rajendran N, Balasubramanian M, Amunts K, Zilles K, Schwartz EL, Fischl B, Triantafyllou C. Locating the functional and anatomical boundaries of human primary visual cortex. NeuroImage. 2009;46(4):915–22. [PMC free article] [PubMed]
Hinds O, Polimeni JR, Rajendran N, Balasubramanian M, Wald LL, Augustinack JC, Wiggins G, Rosas HD, Fischl B, Schwartz EL. The intrinsic shape of human and macaque primary visual cortex. Cereb Cortex. 2008;18(11):2586–95. [PMC free article] [PubMed]
Horner JL, Gianino PD. Phase-only matched filtering. Appl Opt. 1984;23(6):812–6. [PubMed]
Horton JC, Hedley-Whyte ET. Mapping of cytochrome oxidase patches and ocular dominance columns in human visual cortex. Philos Trans R Soc Lond B Biol Sci. 1984;304(1119):255–72. [PubMed]
Jenkinson M, Bannister P, Brady M, Smith S. Improved optimization for the robust and accurate linear registration and motion correction of brain images. NeuroImage. 2002;17(2):825–41. [PubMed]
Jin T, Kim SG. Cortical layer-dependent dynamic blood oxygenation, cerebral blood flow and cerebral blood volume responses during visual stimulation. NeuroImage. 2008;43(1):1–9. [PMC free article] [PubMed]
Jones SE, Buchbinder BR, Aharon I. Three-dimensional mapping of cortical thickness using Laplace's equation. Hum Brain Mapp. 2000;11(1):12–32. [PubMed]
Keilholz SD, Silva AC, Duyn J, Koretsky AP. The contribution of T2* to spin-echo EPI: implications for high-field fMRI studies. Proc ISMRM. 2005;13:32.
Koopmans PJ, Orzada S, Barth M, Norris DG. Distinguishing pial and laminar gradient-echo BOLD signals at 7 Tesla [Abstract] Proc ISMRM. 2009;17:1558.
Lauwers F, Cassot F, Lauwers-Cances V, Puwanarajah P, Duvernoy H. Morphometry of the human cerebral cortex microcirculation: general characteristics and space-related profiles. NeuroImage. 2008;39(3):936–48. [PubMed]
Leite FP, Vanduffel W, Tootell RB, Mandeville JB. Point-image fMRI experiments in awake, behaving macaques. Proc ISMRM. 2004;12:1079.
LeVay S, Hubel DH, Wiesel TN. The pattern of ocular dominance columns in macaque visual cortex revealed by a reduced silver stain. J Comp Neurol. 1975;159(4):559–76. [PubMed]
Logothetis NK. Can current fMRI techniques reveal the micro-architecture of cortex? Nature. 2000;3(5):413. [PubMed]
Menon RS, Goodyear BG. Submillimeter functional localization in human striate cortex using BOLD contrast at 4 Tesla: implications for the vascular point-spread function. Magn Reson Med. 1999;41(2):230–5. [PubMed]
Merboldt KD, Finsterbusch J, Frahm J. Reducing inhomogeneity artifacts in functional MRI of human brain activation-thin sections vs gradient compensation. J Magn Reson. 2000;145(2):184–91. [PubMed]
Moon CH, Fukuda M, Park SH, Kim SG. Neural interpretation of blood oxygenation level-dependent fMRI maps at submillimeter columnar resolution. J Neurosci. 2007;27(26):6892–902. [PubMed]
Ogawa S, Menon RS, Tank DW, Kim SG, Merkle H, Ellermann JM, Ugurbil K. Functional brain mapping by blood oxygenation level-dependent contrast magnetic resonance imaging. A comparison of signal characteristics with a biophysical model. Biophys J. 1993;64(3):803–12. [PubMed]
Olman CA, Inati S, Heeger DJ. The effect of large veins on spatial localization with GE BOLD at 3 T: Displacement, not blurring. NeuroImage. 2007;34(3):1126–35. [PubMed]
Polimeni JR, Balasubramanian M, Schwartz EL. Multi-area visuotopic map complexes in macaque striate and extra-striate cortex. Vision Res. 2006;46(20):3336–3359. [PMC free article] [PubMed]
Polimeni JR, Fischl B, Greve DN, Wald LL. Laminar analysis of high isotropic resolution BOLD activation with a resolution pattern stimulus in human V1 at 7T [Abstract] Proc ISMRM. 2009a;17:1559.
Polimeni JR, Fischl B, Greve DN, Wald LL. Laminar analysis of high isotropic resolution BOLD activation with a resolution pattern stimulus in human V1 at 7T [Abstract] NeuroImage. 2009b;47(Suppl 1):S196.
Polimeni JR, Hinds OP, Balasubramanian M, van der Kouwe A, Wald LL, Dale AM, Fischl B, Schwartz EL. The human V1-V2-V3 visuotopic map complex measured via fMRI at 3 and 7 Tesla [Abstract] Soc Neurosci 2005
Polimeni JR, Hinds O, Triantafyllou C. Spatial scale tuning maps in human visual cortex [Abstract] NeuroImage. 2008;41(Suppl 1):709.
Poser BA, Norris DG. Fast spin echo sequences for BOLD functional MRI. MAGMA. 2007;20(1):11–7. [PMC free article] [PubMed]
Putnam NM, Hofer HJ, Doble N, Chen L, Carroll J, Williams DR. The locus of fixation and the foveal cone mosaic. J Vis. 2005;5(7):3, 632–9. [PubMed]
Ress D, Glover GH, Liu J, Wandell B. Laminar profiles of functional activity in the human brain. NeuroImage. 2007;34(1):74–84. [PubMed]
Shmuel A, Yacoub E, Chaimow D, Logothetis NK, Ugurbil K. Spatio-temporal point-spread function of fMRI signal in human gray matter at 7 Tesla. Neuroimage. 2007;35(2):539–52. [PMC free article] [PubMed]
Schwartz EL. Computational anatomy and functional architecture of striate cortex: a spatial mapping approach to perceptual coding. Vision Res. 1980;20(8):645–69. [PubMed]
Ségonne F, Dale AM, Busa E, Glessner M, Salat D, Hahn HK, Fischl B. A hybrid approach to the skull stripping problem in MRI. NeuroImage. 2004;22(3):1060–75. [PubMed]
Ségonne F, Grimson E, Fischl B. A genetic algorithm for the topology correction of cortical surfaces. Inf Process Med Imaging. 2005;19:393–405. [PubMed]
Silva AC, Koretsky AP. Laminar specificity of functional MRI onset times during somatosensory stimulation in rat. Proc Natl Acad Sci U S A. 2002;99(23):15182–7. [PubMed]
Smirnakis SM, Schmid MC, Weber B, Tolias AS, Augath M, Logothetis NK. Spatial specificity of BOLD versus cerebral blood volume fMRI for mapping cortical organization. J Cereb Blood Flow Metab. 2007;27(6):1248–61. [PubMed]
Tootell RB, Silverman MS, Switkes E, Valois DRL. Deoxyglucose analysis of retinotopic organization in primate striate cortex. Science. 1982;218(4575):902–4. [PubMed]
Triantafyllou C, Hoge RD, Krueger G, Wiggins CJ, Potthast A, Wiggins GC, Wald LL. Comparison of physiological noise at 1.5 T, 3 T and 7 T and optimization of fMRI acquisition parameters. NeuroImage. 2005;26(1):243–50. [PubMed]
van der Kouwe A, Benner T, Salat DH, Fischl B. Brain morphometry with multiecho MPRAGE. NeuroImage. 2008;40(2):559–69. [PMC free article] [PubMed]
Van Essen DC, Maunsell JH. Two-dimensional maps of the cerebral cortex. J Comp Neurol. 1980;191(2):255–81. [PubMed]
Weber B, Keller AL, Reichold J, Logothetis NK. The microvascular system of the striate and extrastriate visual cortex of the macaque. Cereb Cortex. 2008;18(10):2318–30. [PubMed]
Welker W. Why does cerebral cortex fissure and fold? A review of determinants of gyri and sulci. In: Jones EG, Peters A, editors. Cerebral Cortex. 8b. Plenum Press; New York: 1990. pp. 3–136. chapter 10.
Wiggins GC, Wiggins CJ, Potthast A, Alagappan V, Kraft O, Reykowski A, Wald LL. A 32 channel receive-only head coil and detunable transmit birdcage coil for 7 Tesla brain imaging [Abstract] Proc ISMRM. 2006;14:415.
Woolrich MW, Ripley BD, Brady M, Smith SM. Temporal autocorrelation in univariate linear modeling of FMRI data. NeuroImage. 2001;14(6):1370–86. [PubMed]
Yacoub E, Shmuel A, Pfeuffer J, Van De Moortele PF, Adriany G, Andersen P, Vaughan JT, Merkle H, Ugurbil K, Hu X. Imaging brain function in humans at 7 Tesla. Magn Reson Med. 2001a;45(4):588–94. [PubMed]
Yacoub E, Shmuel A, Pfeuffer J, Van De Moortele PF, Adriany G, Ugurbil K, Hu X. Investigation of the initial dip in fMRI at 7 Tesla. NMR Biomed. 2001b;14(7–8):408–12. [PubMed]
Yacoub E, Duong TQ, Van De Moortele PF, Lindquist M, Adriany G, Kim SG, Ugurbil K, Hu X. Spin-echo fMRI in humans using high spatial resolutions and high magnetic fields. Magn Reson Med. 2003;49(4):655–64. [PubMed]
Yacoub E, Shmuel A, Logothetis N, Ugurbil K. Robust detection of ocular dominance columns in humans using Hahn Spin Echo BOLD functional MRI at 7 Tesla. NeuroImage. 2007;37(4):1161–77. [PMC free article] [PubMed]
Yacoub E, Ugurbil K, Olman CA. Feasibility of detecting differential layer specific activations in humans using SE BOLD fMRI at 7 T [Abstract] Proc ISMRM. 2009;17:1607.
Zhao F, Wang P, Hendrich K, Ugurbil K, Kim SG. Cortical layer-dependent BOLD and CBV responses measured by spin-echo and gradient-echo fMRI: insights into hemodynamic regulation. NeuroImage. 2006;30(4):1149–60. [PubMed]
Zheng D, LaMantia AS, Purves D. Specialized vascularization of the primate visual cortex. J Neurosci. 1991;11(8):2622–9. [PubMed]