Abstract
Whether conventional gradient-echo (GE) blood oxygenation-level-dependent (BOLD) functional magnetic resonance imaging (fMRI) is able to map submillimeter-scale functional columns remains debatable mainly because of the spatially nonspecific large vessel contribution, poor sensitivity and reproducibility, and lack of independent evaluation. Furthermore, if the results from optical imaging of intrinsic signals are directly applicable, regions with the highest BOLD signals may indicate neurally inactive domains rather than active columns when multiple columns are activated. To examine these issues, we performed BOLD fMRI at a magnetic field of 9.4 tesla to map orientation-selective columns of isoflurane-anesthetized cats. We could not convincingly map orientation columns using conventional block-design stimulation and differential analysis method because of large fluctuations of signals. However, we successfully obtained GE BOLD iso-orientation maps with high reproducibility (r = 0.74) using temporally encoded continuous cyclic orientation stimulation with Fourier data analysis, which reduces orientation-nonselective signals such as draining artifacts and is less sensitive to signal fluctuations. We further reduced large vessel contribution using the improved spin-echo (SE) BOLD method but with overall decreased sensitivity. Both GE and SE BOLD iso-orientation maps excluding large pial vascular regions were significantly correlated to maps with a known neural interpretation, which were obtained in contrast agent-aided cerebral blood volume fMRI and total hemoglobin-based optical imaging of intrinsic signals at a hemoglobin iso-sbestic point (570 nm). These results suggest that, unlike the expectation from deoxyhemoglobin-based optical imaging studies, the highest BOLD signals are localized to the sites of increased neural activity when column-nonselective signals are suppressed.
Introduction
Conventional gradient-echo (GE) blood oxygenation-level-dependent (BOLD) functional magnetic resonance imaging (fMRI) has become the dominant methodology for brain studies at supramillimeter resolution. However, in-depth investigations of cortical information processing require submillimeter-resolution fMRI to map functional modules such as orientation columns (Hubel and Wiesel, 1962; Blasdel and Salama, 1986). Iso-orientation maps have been obtained by fMRI techniques based on the metabolically related early negative BOLD dip (Kim et al., 2000), hemodynamic cerebral blood flow (CBF) (Duong et al., 2001), or cerebral blood volume (CBV) (Zhao et al., 2005) but not with GE BOLD fMRI (Duong et al., 2000b; Kim et al., 2000). Because the BOLD dip and CBF-based fMRI have poor sensitivity and because the CBV-based fMRI uses exogenous contrast agents, they cannot be easily used for high-resolution human fMRI studies. Additionally, because the commonly used GE BOLD technique is highly sensitive to less-specific large venous vessels (Kim et al., 2000), suppressing draining vessel contributions to positive BOLD (i.e., hyperoxygenated) signal is essential for detecting small neurally activated signals.
Although the draining artifacts can be suppressed and BOLD columnar maps can be obtained (Menon et al., 1997; Menon and Goodyear, 1999; Dechent and Frahm, 2000; Cheng et al., 2001; Goodyear and Menon, 2001), interpreting these maps requires evidence that the BOLD signals are spatially correlated with sites of increased neural activity (Duong et al., 2000b; Cheng et al., 2001). The BOLD signal point spread function (PSF) is determined by a mismatch between the cerebral metabolic rate of oxygen (CMRO2) and CBF responses (Fig. 1). When multiple columns are active, the highest BOLD signals could originate from neurally active domains if CBF PSF is narrow relative to an intercolumn distance (Duong et al., 2001), or the highest BOLD signal could originate from neurally inactive domains if the CBF response is broad as suggested by optical spectroscopic imaging (Malonek and Grinvald, 1996). Therefore, it is critical to determine the neural interpretation of positive BOLD fMRI maps at submillimeter columnar resolution.
In the present study, we examined active column-specific and nonspecific BOLD signals using the cat iso-orientation columns. We also aimed to improve the sensitivity of orientation-selective BOLD signals and to interpret the BOLD iso-orientation maps. To reduce large orientation-nonspecific signals, we evaluated two approaches on continuously acquired data at 9.4 tesla: a continuous stimulation paradigm and spin-echo (SE) BOLD. The continuous cyclic stimulation paradigm with the Fourier analysis approach (Engel et al., 1994; Sereno et al., 1995) is known to be effective for suppressing large orientation-nonspecific signals to detect small orientation-specific signals (Kalatsky and Stryker, 2003); we compared this continuous temporally encoded method to conventional block-design method with differential data analysis. The SE BOLD technique at high magnetic fields, which is known to be mostly sensitive to microvessels (Lee et al., 1999; Zhao et al., 2004), was used for comparison with GE BOLD fMRI. For neural interpretation, BOLD iso-orientation maps were directly compared with CBV fMRI maps and to CBV-weighted optical intrinsic signal (OIS) maps, where spatial neural correlations are known (Fukuda et al., 2006a).
Materials and Methods
Animal preparation.
All experiments were approved by the Institutional Animal Care and Use Committee of the University of Pittsburgh. Six cats (1.3–2.5 kg; 13.6–21.7 weeks of age) were used for fMRI studies, and three of these cats were studied by optical imaging for comparison. Each cat was initially treated with atropine sulfate (0.05 mg/kg, i.m.). Anesthesia was induced by a mixture of ketamine (10–20 mg/kg, i.m.) and xylazine (1 mg/kg, i.m.); surgery was performed under a mixture of air and O2 (maintaining a 30–35% O2 level) with 2% isoflurane, and all experiments were performed at a 0.8–1.2% isoflurane level. The cat was paralyzed with pancuronium bromide (0.2 mg/kg/hr, i.v.), and appropriate contact lenses were used. For fMRI studies, the cat was placed in a cradle and secured in a normal postural position by a custom-designed head frame. After fMRI experiments, the cat was moved to the room adjacent to the MR system for the purpose of optical imaging and secured in a stereotaxic apparatus. A cranial window was made over the location of the MRI coil, and the exposed cortical surface was covered by agarose (∼1% in saline). Throughout all experiments, rectal temperature and end-tidal CO2 were maintained at 37.9 ± 0.2°C and 3.5 ± 0.1%, respectively, as in typical anesthetized cat experiments but slightly lower than the awake condition (Herbert and Mitchell, 1971).
Visual stimulation.
To activate multiple iso-orientation columns (Fig. 1), high-contrast square-wave full-field moving gratings (0.15 cycle/°; 2 cycles/s; moving direction reversal per 0.5 s) (Movshon et al., 1978) were presented binocularly. Visual stimulation was generated using either software from Cogent Graphics (John Romaya, The Wellcome Department of Imaging Neuroscience, University College London, London, UK; http://www.vislab.ucl.ac.uk/Cogent/index.html) or VSG2/5 (Cambridge Research Systems, Rochester, UK). For each continuous stimulation run, one stimulation cycle consisting of eight consecutive orientations (0–157.5°; 22.5° increments; 10 s each) was repeated 11 times without gaps between stimuli (80 s per cycle; 880 s per run). To obtain functional responses in a single-orientation stimulation paradigm, 20 s one-orientation stimulation (67.5°) was presented at 50 s before and 110 s after each continuous stimulation run. For comparison, a block-design stimulation run consisting of four 10 s orientations (0, 45, 90, and 135°) alternated with 20 s homogeneous gray (“blank” control) was performed (10 cycles; 1200 s) using GE BOLD fMRI. Stimulation was synchronized with the image acquisition.
MRI experiments.
All MR imaging was performed on a 9.4 tesla MR system (Varian, Palo Alto, CA) with a horizontal magnet (clear bore size, 31 cm) and an actively shielded gradient coil (inner diameter, 12 cm; maximal strength, 40 gauss/cm; rise time, 130 μs) using a custom-made surface radio frequency coil (inner diameter, 1.6 cm) positioned over the primary visual cortex. The position of the functional imaging slice was determined based on a flow-compensated, gradient-recalled, three-dimensional (3-D) venographic image [repetition time (TR), 50 ms; echo time (TE), 20 ms; data matrix, 512 × 256 × 256; field of view (FOV), 3.5 × 2.2 × 2.2 cm3; resolution, 68 × 86 × 86 μm3] (Park and Kim, 2005). Using the 3-D venographic image, a single 1-mm-thick imaging slice was selected with the aid of a homemade C++ visualization program (Microsoft, San Francisco, CA) based on two criteria: (1) avoiding large surface veins that induce large susceptibility artifacts in echo planar imaging (EPI) and may have large nonspecific GE BOLD signal changes (Kim et al., 2000); and (2) including a flat dorsal surface area as large as possible and tangential to the marginal gyrus. In this slice, the intracortical veins perpendicular to the surface of gyrus appear as dark spots (Fig. 2 B) (columns are supposed to be arranged parallel to these intracortical veins). Once a slice position was determined, a two-dimensional (2-D) vessel-weighted GE image was acquired (TR, 50 ms; TE, 15–20 ms; data matrix, 256 × 256; FOV, 2 × 2 cm2; thickness, 1 mm) as an anatomical reference of fMRI data.
Three different fMRI images were obtained using EPI techniques with data matrix, 64 × 64, and FOV, 2 × 2 cm2 (thickness, 1 mm): (1) GE BOLD fMRI with TR = 0.5 s and TE = 18 ms; (2) SE BOLD fMRI with TR = 2.0 s and TE = 40 ms (Lee et al., 1999); and (3) CBV-weighted fMRI with TR = 1.0 s and TE = 10 ms after an intravascular bolus injection of a dextran-coated monocrystalline iron oxide nanoparticles (MION) contrast agent (10–20 mg Fe/kg body weight) (Zhao et al., 2005, 2006; Fukuda et al., 2006a). Average repetitions of continuous stimulation runs (n = 6 cats) were 3.3 ± 1.8 for GE BOLD, 5.7 ± 1.0 for SE BOLD, and 2.5 ± 0.5 for CBV-weighted fMRI.
OIS experiments.
After the fMRI experiments, all OIS imaging with 570 ± 10 nm wavelength light (n = 3 cats; three runs) was performed on a custom-built imaging system (Moon et al., 2004). The 570 nm wavelength was chosen to record CBV-weighted OIS because of its high sensitivity compared with deoxyhemoglobin-weighted 620 nm OIS when continuous stimulation is used. Temporally encoded orientation maps obtained from these two wavelengths were almost identical [Fukuda et al. (2006a), their supplemental Fig. 2]. The camera (FOV, 1.9 × 1.4 cm2; 29 × 29 μm2 per pixel) was positioned above the MR imaging area and adjusted so that the pial vessel pattern closely matched that of the MR image. The focus plane was 200–500 μm below the cortical surface. Images were acquired with a temporal resolution of 33.3 ms.
Data analysis.
All data were analyzed with homemade programs coded by Visual C++ (Microsoft, Seattle, WA) or Matlab (The MathWorks, Natick, MA). Quantitative comparisons were performed on 10 hemispheres in six cats for fMRI studies and four hemispheres in three cats for fMRI versus OIS studies.
Preprocessing.
All MR images were interpolated by zero-padding in k-space from 64 × 64 matrix to 128 × 128 matrix, resulting in nominal in-plane resolution of 156 × 156 μm2 per pixel. Before zero-filling, a low-pass Hanning filter with a cutoff frequency of 1.60 mm−1 (Oppenheim and Schafer, 1989) was applied to avoid the ringing artifact that resulted from subsequent zero-filling. This procedure did not change activation patterns (supplemental Fig. 1, available at www.jneurosci.org as supplemental material).
Because the initial 80 s (one-cycle data after the onset of continuous stimulation) included the time to reach steady state, the following continuous stimulation (800 s, 10-cycle) data were selected for additional analysis. To remove slow signal fluctuations observed in all brain pixels (referred to as global signal), signal intensity in each pixel at each time point was divided by the global signal intensity at the corresponding time point. The global signal was calculated by averaging over all brain pixels (with intensity higher than 80% of the mean intensity of the image). The underlying assumption of the global normalization is that the orientation-specific signal within the large region of interest (ROI) remains constant over varying orientation stimulation. Because each pixel has an orientation-specific temporal modulation with a phase corresponding to orientation preference and the phases are distributed evenly from 0 to 2π in the large ROI, the sum of orientation-specific signals of all pixels in the large ROI is close to constant. Therefore, the global normalization procedure did not affect orientation-specific BOLD signals. It should be noted that global signal normalization was applied only to continuous stimulation data, not to block-design stimulation data.
For CBV-weighted fMRI data, we first extracted CBV responses from the CBV-weighted and BOLD fMRI signals by adjusting different echo times (Zhao et al., 2006). Stimulation simultaneously induces CBV increase and deoxyhemoglobin (dHb) content decrease (i.e., hyperoxygenation); an increase in MION content induced by CBV increase decreases fMRI signals, whereas a decrease in dHb content increases fMRI signals (i.e., positive BOLD response). To determine only stimulation-induced CBV responses without the contribution of dHb change, BOLD contribution to the CBV-weighted fMRI signals was corrected as (CBV-weighted response) minus (BOLD response) × (10 ms/18 ms), where 10 ms and 18 ms are TE of CBV-weighted and BOLD measurements, respectively. This correction assumes that the intravascular contribution to BOLD signals is minimal in our experimental conditions. Columnar patterns of fMRI signals with MION injection did not change before and after BOLD correction (r = 0.96 ± 0.01; n = 10 hemispheres; p < 0.001 for each) (supplemental Fig. 2, available at www.jneurosci.org as supplemental material), because the effect of MION in vessels almost overcomes the BOLD effect.
Generation of functional maps.
For block-design stimulation data, single-condition maps were obtained by subtracting an average of 5 s prestimulus baseline images from an average of images obtained between 5 and 10 s after the onset of each orientation stimulus (i.e., 0, 45, 90, and 135°). Differential maps were determined by subtracting the single-condition maps obtained from two orthogonal 0 and 90° stimuli.
For continuous stimulation data, iso-orientation maps were determined by two different approaches: differential method and Fourier analysis (Engel et al., 1997; Kalatsky and Stryker, 2003). For the differential method, an average image obtained at 5–10 s after the onset of 90° stimulation was subtracted from that corresponding to 0° stimulation. To generate temporally encoded iso-orientation maps with Fourier analysis, the orientation-specific fMRI signal induced by continuous stimulation was modeled using a sinusoidal function with magnitude and phase at an orientation-specific frequency, f s, as follows: where ΔS(x, y, t) is the orientation-specific signal intensity at position x and y at time t (i.e., “activation” map or iso-orientation map); M(x, y) and φ(x, y) are the absolute magnitude and phase of the frequency f s component [i.e., 1/(80 s) = 0.0125 Hz for orientation-specific stimulation cycle]. Phase φ(x, y) is a sum of the stimulation onset time and a slow hemodynamic delay. Therefore, the hemodynamic delay must be considered for assignment of orientation preference to each pixel; it was found to be ∼5 s for BOLD fMRI and ∼13 s for CBV fMRI and OIS data (Fukuda et al., 2006a) (supplemental Fig. 3, available at www.jneurosci.org as supplemental material). Thus, the iso-orientation map of 0° orientation presented at t of 0–10 s was obtained at t of 5 s for BOLD and t of 13 s for CBV studies. Because the orientation-specific signal intensity is dependent on the sensitivity of the method used, it is normalized by corresponding mean magnitude (M̄) of M(x, y) within the region of interest for comparisons between different activation maps. Note that iso-orientation maps obtained from Fourier analysis are similar to “differential” maps.
Stimulation-induced hyperoxygenation increases the BOLD signal from its baseline because of decreases of dHb amount, whereas an increase in CBV decreases CBV-weighted fMRI signal and OIS from their baselines. Thus, the polarity of the BOLD fMRI signal is opposite to both CBV-weighted fMRI signal and OIS. To match the polarity of BOLD and CBV maps, values of CBV-weighted fMRI and OIS maps were inverted. Thus, brighter pixels in functional maps indicate higher positive BOLD signals or higher CBV changes. It should be noted that the polarity of the CBV signal is not inverted when its time course is presented in figures.
All iso-orientation maps were spatially filtered using a bandpass filter with a cutoff-frequency of 0.20–1.66 mm−1 to remove high-frequency noise and low spatial frequency change (supplemental Fig. 1, available at www.jneurosci.org as supplemental material). To coregister between MRI and OIS images, commonly identifiable vessels were selected in anatomic MR and optical images, and then their positions were coregistered by the linear affine transformation (Zwillinger, 1995). The same transformation was applied to functional maps.
Quantitative analyses.
All quantitative analyses [percentage signal change, reproducibility, intercolumn distance, contrast-noise ratio (CNR), spatial specificity ratio, and similarity] were performed on an active ROI (17–33 mm2) in each hemisphere, which was determined based on an activation map of CBV-weighted fMRI or OIS. As a control, similar analyses were also performed on a noise ROI (20–35 mm2), which consists of a similar number of pixels used for the active ROI but around the four corners of the image.
To obtain percentage changes for block-design stimulation data, the signal intensity in single-condition maps was normalized with the average prestimulus baseline intensity. For the continuous stimulation data, the modeled signal with a sinusoidal function (see Eq. 1) oscillates between −M(x, y) and M(x, y), thus an orientation-specific response is 2 × M(x, y). Then, percentage changes were calculated from dividing 2 × M(x, y) of the orientation-specific frequency by the magnitude of zero frequency component.
Reproducibility was examined cycle-by-cycle and run-by-run. For cycle-by-cycle reproducibility, either differential or 0° iso-orientation maps for each cycle were compared with the averaged 0° iso-orientation map obtained from all cycles using linear correlation analysis. For run-by-run reproducibility, an iso-orientation map was obtained from each 800 s run using Fourier analysis, and then a pixel-wise correlation between the maps obtained from the first run and subsequent runs was calculated and averaged for all eight orientations.
To calculate an average distance between iso-orientation columns, activation foci were determined, and then the average distance of neighboring foci was determined by the nearest-neighbor spatial analysis technique (Murphy et al., 1998) within the range of interpeak distances between 313 μm and 2 mm. The CNR of the orientation-selective signal was calculated by dividing the signal magnitude at 0.0125 Hz in the active ROI by that in the noise ROI.
Hemodynamic response induced by orientation-selective stimulation in the primary visual cortex consists of small orientation-specific and large orientation-nonspecific components. The ratio of orientation-specific to total (total = orientation-specific + orientation-nonspecific) response can be used to evaluate spatial specificity of the signals [Fukuda et al. (2006b), their Fig. 8]. Briefly, the ratio will be close to 1.0 if the orientation-specific signal is dominant in the total signal (i.e., narrow hemodynamic PSF in Fig. 1), whereas the ratio will be close to 0.0 if the orientation-nonspecific signal is dominant in the total signal (i.e., broad PSF in Fig. 1). For continuous stimulation data, the orientation-nonspecific response is saturated because of continuous cyclic stimulation without gaps and cannot be accurately determined. Therefore, the total response induced by orientation-selective stimulation was not obtained from continuous stimulation data but instead from the 20 s one-orientation stimulation data before the initiation of continuous stimulation.
To determine the similarity of iso-orientation maps obtained from different fMRI methods, normalized iso-orientation maps [ΔS(x, y, t)/M̄] of the same stimulus orientation were compared using a linear correlation analysis, and correlation coefficient values (R) for all eight stimulus orientations were averaged. To quantify the similarity of columnar patterns between BOLD fMRI and OIS iso-orientation maps, the least-distance for local maxima of iso-orientation domains between two 0° iso-orientation maps was determined.
Statistical analyses.
All data were analyzed by Student's t test (paired, two-tailed) or one-way ANOVA. The p value of the correlation coefficient was determined after Bonferroni correction. Data are presented as mean ± 1 SD.
Results
Orientation-specific signal of conventional BOLD fMRI
To reconfirm that orientation columns are very difficult to be detected with standard block-design stimulation (Duong et al., 2000b; Kim et al., 2000), conventional GE BOLD images were acquired during the four-orientation stimulation protocol in a single cat. If the BOLD fMRI signal is highly orientation selective, the activation patterns of orthogonal stimuli should be complementary to each other (Blasdel and Salama, 1986; Blasdel, 1992). In single-condition maps (one orientation vs blank control), activation with an average signal change of 2.4% within the ROIs (Fig. 2 B, black dashed rectangles) was diffused across the whole brain imaging area and similarly inhomogeneous regardless of stimulus orientations (Fig. 2 C, first and second panel) (45 and 135° not shown), indicating that the orientation-nonspecific signal is dominant over the orientation-specific signal. The sites of the largest BOLD activation were located at the midline and at the edges of the brain where large draining vessels exist (Fig. 2 C, black arrowheads). To remove orientation-nonspecific signals, which could mask small orientation-specific modulations, differential maps were obtained from two single-condition maps of orthogonal stimuli (0 vs 90°) (Fig. 2 C, third panel). Signal intensities in the differential map are at a background noise level except for a few patchy areas that correspond to large veins (Fig. 2 C, white arrowheads). Consequently, the cycle-by-cycle reproducibility of differential maps was quite poor in the active ROI, similar to noise (Fig. 2 D). To further explore signal properties, regions with higher signals responding to 0 or 90° stimulus were determined from the differential map, and time courses within 0 and 90° ROIs (Fig. 2 E, first panel) were obtained during 0 and 90° stimulation (Fig. 2 E, second and third panels). Clearly, the difference of signal changes during 0 and 90° stimulation (green traces) was very small relative to total stimulation-induced responses (red and blue traces). A cycle-by-cycle variation of single-condition signals for 5–10 s after stimulation onset was ∼2.5–4.5% (see error bars in red and blue traces), whereas the average intensity of the differential signal (green traces) was ∼0.3–0.7%. The average differential signal was also smaller than variations of the differential signals: 0.28–0.44 times SDs of differential signals across 10 cycles. As a result, large variations of both single-condition and differential signals will reduce the detectability of orientation-selective signals. It should be noted that unlike BOLD signals in the previous study (Kim et al., 2000), the early dip was not observed in the present study. One obvious discrepancy between the previous and present studies is temporal dynamics of positive BOLD signals [time-to-peak, 8–10 s (Kim et al., 2000) vs ∼5 s (in present studies)]. The magnitude of the early dip is closely dependent on the animal's condition [Harel et al. (2002), their Fig. 4]; when a relatively higher end-tidal CO2 level was used, stimulation-induced hemodynamic response is delayed (Cohen et al., 2002), resulting in a higher probability to detect the early dip. Note that the differential method robustly detected the ocular dominance column-specific signals in humans, the amplitude of which seems larger, ∼0.6–1.2% [TE, 15 ms at 4T (Cheng et al., 2001)], than our orientation-specific signal, 0.3–0.7% (TE, 18 ms at 9.4 T).
To effectively reduce dominant, unstable orientation-nonspecific signals, continuous stimulation without blank control was adopted. Global signal fluctuations over a whole brain imaging area (see above, Data analysis) is likely because of changes in animal physiology, including blood pressure and blood gases (Fig. 3 A, black trace). Similar signal fluctuation observed in the small ROI (red trace) masks a small orientation-specific signal (i.e., a modulation of 80 s cycle). Thus, signal intensity in each pixel at each time point was divided by the global signal intensity at the corresponding time point to remove the slow fluctuation of baseline signal over a long experiment (blue trace). To examine the effects of continuous stimulation, cycle-by-cycle reproducibility of differential maps was calculated within the same ROI shown in Figure 2 B. The continuous stimulation paradigm dramatically improved the sensitivity of orientation-specific signals compared with the block-design stimulation paradigm (compare Figs. 3 B and 2 D); average correlation coefficients across 10 cycles were 0.44 ± 0.17 for continuous stimulation data versus 0.28 ± 0.29 for block-design stimulation data.
To explore different methods of analyzing the same data obtained with continuous temporally encoded stimulation, normalized data of the 80 s stimulus cycle were averaged over 10 cycles and signal intensities from the small ROI were plotted as a function of time (Fig. 3 C, blue dots). The modulation induced by a change of stimulation orientation is well approximated by a sinusoidal function (Fig. 3 C, brown trace) because of the smooth change in orientation selectivity of neurons in primary visual cortex (Swindale, 1998) and the linear transform characteristics of hemodynamic responses (Boynton et al., 1996) and can easily be extracted by a frequency analysis method (Engel et al., 1997; Kalatsky and Stryker, 2003). The cycle-by-cycle reproducibility of orientation-selective maps from 800 s continuous stimulation data analyzed with Fourier analysis at the orientation-selective frequency (1/80 s−1) is highly significant within the ROI shown in Figure 2 B (Fig. 3 D). This demonstrates that the Fourier analysis of continuous temporally encoded stimulation data further improves the sensitivity of orientation-specific signals; average correlation coefficients across 10 cycles were 0.57 ± 0.11 with Fourier analysis (supplemental Fig. 4, available at www.jneurosci.org as supplemental material). The spatial filter did not change the shapes of the isoorientation patches (supplemental Fig. 1, available at www.jneurosci.org as supplemental material), but it did improve the spatial homogeneity of the orientation map, which resulted in better reproducibility; the average correlation coefficient for eight orientation stimulations across 10 runs with filter (0.60 ± 0.07) was significantly higher (p = 6.66 × 10−43) than without filter (0.57 ± 0.11). Thus, the continuous temporally encoded stimulation paradigm and Fourier analysis with a spatial filter can be used to obtain temporally encoded iso-orientation maps.
Temporally encoded BOLD orientation columnar mapping
To test the reproducibility of temporally encoded iso-orientation maps obtained using conventional GE BOLD fMRI across different runs, the maps of 0° stimulation determined from three consecutive 800 s runs were compared (Fig. 4 A). Bright pixels indicate increased BOLD signals during 0° stimulation, whereas dark pixels correspond to increased BOLD signals during 90° stimulation (i.e., decreased signal during 0° stimulation relative to the mean intensity). The average interdistance of iso-orientation patches for all eight orientations was 1.40 ± 0.05 mm (n = 10 hemispheres), which is consistent with previous reports of 1.1–1.4 mm (Lowel and Singer, 1990; Duong et al., 2000b, 2001; Kim et al., 2000; Zhao et al., 2005; Fukuda et al., 2006a). Patchy patterns were quite reproducible across three 800 s runs (Fig 4 A, yellow or black dashed regions); pixel-wise correlation coefficients (Fig. 4 B) within the left hemisphere ROI (Fig. 4 A, yellow dashed contour) are significantly high (r = 0.78 at first vs second runs and 0.87 at first vs third runs; p < 0.001 for each). The average pixel-wise correlation coefficient for all eight orientations across three runs was 0.74 ± 0.12 (p < 0.001 for each hemisphere) for 10 hemispheres, which is similar to the reproducibility of iso-orientation maps obtained by CBV-weighted fMRI (Fukuda et al., 2006a) (r = 0.75 ± 0.10; n = 5 at first vs second runs with 128 × 128 data collection and 512 × 512 interpolation). It should be noted that the correlation coefficients in noise regions were close to zero (r = −0.01 ± 0.06 at first vs second runs and 0.00 ± 0.12 at first vs third runs; n = 6 cats).
To further examine the reproducibility of BOLD maps, iso-orientation maps were compared with the maps obtained 1 week later (Fig. 5 A) and 1 month later (Fig. 5 B). Because imaging slice positions were varied at different sessions, images were coregistered based on intracortical vessels in anatomical images (Fig. 5, blue rectangles). Most patches appeared in the same locations even 1 week and 1 month later (Fig. 5, plus signs). The results suggest that the temporally encoded BOLD iso-orientation map is highly reproducible.
Comparison of GE and SE BOLD orientation-specific activation
Despite the high reproducibility of iso-orientation maps obtained from conventional BOLD signals, the maps may be distorted by stimulation-related draining signals, because the GE BOLD signal is extremely sensitive to large draining veins (Lai et al., 1993; Menon et al., 1993; Frahm et al., 1994; Kim et al., 1994, 2000; Lee et al., 1999; Duong et al., 2000a,b). To examine this issue, SE BOLD fMRI (Lee et al., 1999; Zhao et al., 2004) was performed during exactly the same continuous stimulation paradigm, and the GE and SE BOLD magnitude maps of orientation-specific signals were compared. Artifacts from large draining vessels at the edges and at the midline of the brain were determined by the single-condition GE BOLD map obtained from block-design stimulation (Fig. 6 A,B, white arrowed regions). The magnitude of orientation-specific signals in continuous stimulation data are ∼20% of the signal intensity in single-condition maps because of the reduction of orientation-nonspecific contributions. However, the high signal intensity regions observed in the single-condition map still remained at the edges and midline of the brain in the GE BOLD magnitude map but not in the SE BOLD magnitude map (Fig. 6 C,D, white arrowed regions). Thus, large vessel contributions can distort GE BOLD-based functional maps obtained even from continuous stimulation and Fourier analysis. In the internal cortical regions excluding large pial vessel areas (within the magenta rectangles in Fig. 6 C,D), the orientation-specific activation patterns are similar between the normalized GE and SE BOLD fMRI magnitude maps (Fig. 6 E); a pixel-by-pixel correlation between these two maps is significantly high (r = 0.62; p < 0.001) (Fig. 6 F). A major difference between GE and SE BOLD signals in the internal cortical region was sensitivity. Without normalization by M̄, the actual signal intensity of orientation-specific signals in the internal cortical region for GE BOLD (0.31 ± 0.10%; n = 10 hemispheres) was significantly higher (p < 0.001) than that for SE BOLD fMRI (0.19 ± 0.07%; n = 10). The mean CNR of averaged runs for the orientation-specific GE BOLD signal (7.39 ± 2.66; n = 10; 3.3 runs) was significantly higher (p < 0.001) than that for the orientation-specific SE BOLD signal (3.09 ± 0.99; n = 10; 5.7 runs); CNRs of single runs for GE and SE BOLD signals were 4.73 ± 1.65 and 1.78 ± 0.40 (n = 10), respectively.
To obtain insights into the orientation-nonselective contribution of parenchymal vessels on BOLD signals, the spatial specificity ratio (see Materials and Methods) was obtained. The spatial specificity ratio of GE BOLD fMRI within the internal cortical region (0.17 ± 0.06; n = 10 hemispheres) was significantly lower (p = 0.007) than that of SE BOLD fMRI (0.26 ± 0.09; n = 10), suggesting that GE BOLD signal is more spread than the SE BOLD signal, possibly because of artifacts in and around draining intracortical veins.
Assignment of stimulus orientation to temporally encoded BOLD maps
Higher BOLD signal changes could occur at cortical columns with lower CBF and CBV responses, depending on the relationship between PSFs in oxygen consumption and CBF changes (Fig. 1). To examine this issue, BOLD maps were compared with CBV fMRI maps with the known neural interpretation (Fukuda et al., 2006a). The CBV iso-orientation maps are well matched with GE and SE BOLD maps except for the regions at the edge of the brain (Fig. 7 A, white arrowheads) (for another example, see supplemental Fig. 5, available at www.jneurosci.org as supplemental material). Pixel-wise scatter plots within the active ROI in the right hemisphere (Fig. 7 B) show that both the GE and SE BOLD maps are significantly correlated with the CBV fMRI map (r = 0.83 and 0.86, respectively; p < 0.001 for each), and that the GE BOLD map is significantly correlated with the SE BOLD map (r = 0.82; p < 0.001) (Fig. 7 C). Similar results were obtained from all six cats. Average pixel-wise correlation coefficients within an active ROI across 10 hemispheres are 0.78 ± 0.06 for GE BOLD versus CBV fMRI, 0.76 ± 0.08 for SE BOLD versus CBV fMRI, and 0.72 ± 0.10 for GE versus SE BOLD fMRI maps. One-way ANOVA revealed that the correlation coefficients among iso-orientation maps from the three fMRI techniques have no statistical difference (F (2,29) = 1.386; p = 0.267). These results suggest that, except for the large vessel regions, higher orientation-specific BOLD signals are located at larger orientation-specific CBV response regions and likely mark the sites of increased neural activity.
To further evaluate the assignment of BOLD iso-orientation maps to neuronal activation, optical imaging with a 570 nm wavelength was performed after fMRI in three cats. OIS maps can serve as a reference to a spatial pattern of neuronal activity; a good correlation between neural activity and OIS has already been observed in the visual cortex (Grinvald et al., 1986; Shmuel and Grinvald, 1996; Maldonado et al., 1997; Bosking et al., 2002). Although neural correlates to OIS in visual cortex have been demonstrated only with dHb-weighted wavelengths (e.g., 610–650 nm), good agreement between iso-orientation maps of 570 nm OIS and those of dHb-weighted OIS (Fukuda et al., 2005, 2006a) suggests that 570 nm OIS is also a good neural correlate. To compare BOLD fMRI and OIS orientation maps, vessel patterns in OIS anatomical images were coregistered with those in the MR surface image (Fig. 8 A–C, left panels, green traces). After the coregistration, BOLD fMRI iso-orientation maps generally agree well with the OIS iso-orientation map (Fig. 8 A–C, right panels, plus signs). If the cortical surface of the coregistered region is slightly curved, signals of OIS and fMRI at the coregistered pixels have different anatomical origins. Thus, small differences between BOLD fMRI and OIS orientation maps are expected. The average least distances of local maxima (Fig. 8 A–C, top right panels) for GE BOLD versus OIS (0.22 ± 0.02 mm; n = 4 hemispheres), SE BOLD versus OIS (0.22 ± 0.02 mm), and GE BOLD versus SE BOLD fMRI iso-orientation maps (0.23 ± 0.01 mm) were similar (F (2,11) = 0.673; p = 0.534). Because these distances are shorter than a column width (i.e., a half of the measured intercolumnar distance, ∼0.70 mm), we consider the BOLD iso-orientation maps to match with OIS iso-orientation maps. Thus, the region of highest orientation-specific BOLD signal is likely to indicate the site of increased neural activity.
Discussion
We demonstrated for the first time that highly reproducible iso-orientation maps can be obtained using the conventional positive BOLD fMRI in combination with continuous temporally encoded stimulation and Fourier analysis. Our major finding is that both GE and SE BOLD fMRI orientation maps excluding large pial vessel regions are significantly correlated with CBV fMRI and CBV-weighted OIS orientation maps. Thus, the PSF of BOLD response is sufficiently narrow to resolve iso-orientation columns (Fig. 1, solid red line); the highest orientation-selective BOLD signals indicate the sites of increased neural activity. In other words, the magnitude of BOLD response to preferred stimulation should be larger than the magnitude of response to nonpreferred stimulation (Fig. 9). Because intercolumnar distance of human ocular dominance columns revealed by a positive BOLD signal (Cheng et al., 2001) is similar to the cat iso-orientation columns, our finding suggests that the higher BOLD signal changes induced by left (right) eye stimulation in humans properly indicate left (right) ocular dominance columns.
Mapping orientation columns with BOLD fMRI
We confirmed that a conventional positive BOLD signal responding to single orientation stimulation is widespread and strongly weighted to draining vessels and thus may not be used for single-condition mapping at submillimeter columnar resolution (Kim et al., 2000). Because the orientation-specific signal is much smaller than variations of orientation-nonspecific signals in our study, the differential method cannot effectively separate the orientation-specific signal from the orientation-nonspecific signal. In contrast, the continuous stimulation approach minimizes relatively large time-variable orientation-nonspecific responses during a long data acquisition time, because it almost saturates the orientation-nonspecific signals, including draining artifacts. The continuous stimulation with Fourier data analysis is also expected to increase the detectability or power of the orientation-specific signal as follows. First, the orientation-selective response to the gradual change of stimulus orientation is well fitted by a sinusoidal function (Fig. 3 C) (supplemental Fig. 3, available at www.jneurosci.org as supplemental material). In other words, previous knowledge about a neural tuning curve is needed. If the hemodynamic response is not modulated smoothly over time even with cyclic continuous stimulation, the power of BOLD orientation-specific signal should be significantly reduced (supplemental Fig. 4, available at www.jneurosci.org as supplemental material). Thus, if a single orientation is presented repeatedly as a block-design stimulation paradigm, “single-condition” mapping of orientation columns is extremely difficult with BOLD data collection and Fourier analysis. Second, this orientation-specific stimulation cycle can be set at a frequency far from the physiological periodic components to avoid any contamination of artifacts. Therefore, Fourier analysis in combination with temporally encoded continuous stimulation is more efficient at detecting orientation-selective responses compared with the differential method with the block-design stimulation paradigm (Kim et al., 2000; Duong et al., 2001; Zhao et al., 2005).
Spatial specificity and sensitivity of BOLD techniques
The draining artifact from GE BOLD signal cannot be completely eliminated by continuous stimulation and frequency analysis (Fig. 6), and it is further reduced using the SE BOLD fMRI technique. Consequently, SE BOLD fMRI increases the spatial specificity to active cortical columns (0.26 for SE BOLD vs 0.17 for GE BOLD). Although the draining artifact is negligibly small in the SE BOLD signal, the spatial specificity of the SE BOLD signal is still lower than that of other fMRI signals: an early CMRO2-related signal (0.4–1.0) (Kim et al., 2000), CBF fMRI (0.7) (Duong et al., 2001), and CBV-weighted fMRI (0.4) (Zhao et al., 2005). We believe that the poor spatial specificity of the BOLD signal results from a mismatch of CMRO2 and CBF responses; an increase in CMRO2 response decreases the BOLD signal, whereas an increase in CBF response increases the BOLD signal. Thus, the relative signal difference between active and inactive columns for BOLD signals should be much smaller than that for either CMRO2 or CBF signals.
In addition to specificity, sensitivity is another indispensable factor for utility of brain-mapping techniques. CNR of the GE BOLD signal was approximately three times higher than that of SE BOLD. Thus, to achieve CNR comparable with GE BOLD, the SE BOLD measurement requires a signal averaging ∼32 times more. Similarly, measurements of CBF (Kwong et al., 1992; Edelman et al., 1994; Kim, 1995; Duong et al., 2001) and early negative dip signals (Kim et al., 2000) have poor sensitivity despite their high spatial specificity. Compared with these techniques, the CBV-weighted measurement (Kennan et al., 1998; Mandeville et al., 1998; van Bruggen et al., 1998; Mandeville and Marota, 1999; Kim and Ugurbil, 2003) has the highest sensitivity (CNR 7.15 in CBV vs 4.73 in GE BOLD per run), but this approach cannot be easily used in human studies, because a large amount of exogenous contrast agent is required. Therefore, in practice, the GE BOLD technique would be an appropriate choice for mapping human functional columns within a limited recording time in the parenchymal region excluding the large pial vessels. However, if the recording time is sufficiently long to improve CNR, the SE BOLD technique at high fields should be used because of its reduced sensitivity to large draining vessels.
Orientation-specific BOLD signal source
Our BOLD studies revealed that the preferred stimulation induced a larger positive BOLD signal than the nonpreferred stimulation. Because positive BOLD signal results from a decrease in dHb content resulting from hyperoxygenation, this finding suggests that dHb content decreases more during the preferred stimulation than during the nonpreferred stimulation (i.e., more hyperoxygenation during preferred stimulation). However, this result contradicts the expectation of dHb content changes obtained from optical spectroscopic imaging (Malonek and Grinvald, 1996); the preferred stimulation induced a smaller decrease in the dHb amount than the nonpreferred stimulation (Fig. 9). More hyperoxygenation for nonpreferred stimulation is also consistently observed in OIS at dHb-weighted wavelengths (Grinvald et al., 2000; Fukuda et al., 2005).
In our study, the fMRI signal originates from a 1-mm-thick slice centered on the middle cortical layer, whereas the OIS is from the pia mater and upper cortical layers because of limited light penetration into the cortex (Kennerley et al., 2005). However, this difference of anatomical origins unlikely leads to the opposite polarity of dHb signals between fMRI and optical imaging signals. It should be noted that when the CBF and CBV responses induced by neural activity are minimal at the vasodilator-induced hypotension condition (Nagaoka et al., 2006), both BOLD signal (Fukuda et al., 2006a) and OIS (Fukuda et al., 2006b) have similarly larger hypo-oxygenated changes at active columns than inactive columns, similar to the condition of the early dip. This indicates that the dHb contribution is dominant in both BOLD signal and OIS when CBF and CBV responses are minimal, suggesting that the difference between BOLD signal and OIS at the normal condition appears to relate to CBF and/or CBV responses.
In fMRI studies, the stimulation-induced arterial blood velocity increase results in a large delivery of fresh (unexcited) water into the imaging slice, which gives rise to an increase in the signal change (commonly referred to as “inflow effect”) (Duyn et al., 1994; Frahm et al., 1994; Kim et al., 1994; Gao et al., 1996). In our SE BOLD fMRI, the inflow contribution to fMRI can be <0.1% based on a previous study with the same pulse sequence and similar imaging parameters (Jin et al., 2006). This small amount of inflow contribution could not be ignored because of the extremely small orientation-specific SE BOLD signal (∼0.2%). If the inflow effect significantly contributes to the orientation-specific signal in BOLD measurements, the polarity of the orientation-selective signal may follow that of the CBF response (Duong et al., 2001). However, our preliminary result suggests that the inflow effect is unlikely to change the polarity of the orientation-selective BOLD signal (supplemental Fig. 6, available at www.jneurosci.org as supplemental material).
In OIS studies, when dHb-weighted wavelengths are used, the contribution of oxy-hemoglobin signal to dHb-weighted OIS is not negligible if the CBV response is significant. An increase in CBV makes hypo-oxygenated OIS (i.e., dip) larger and hyperoxygenated OIS smaller because of an increase in light absorption. Additionally, blood velocity changes in microvessels will induce light scattering, which may contribute to the observed OIS data. To separate possible contribution of CBV (total hemoglobin) and light scattering signals from dHb signal, OIS can be decomposed into components of oxy-hemoglobin, deoxy-hemoglobin (total hemoglobin is the sum of oxy- and deoxy-hemoglobin) and light scattering using spectroscopic analysis. Results from decomposed dHb signal were similar to those from dHb-weighted OIS. However, a model used for OIS decomposition in cat and monkey visual cortex (Malonek and Grinvald, 1996; Shtoyerman et al., 2000; Fukuda et al., 2005) was oversimplified (Mayhew et al., 1999; Nemoto et al., 1999; Lindauer et al., 2001; Sato et al., 2002) and may contain a substantial error.
Apparent contradiction between BOLD signal and OIS is also found in the magnitude of early dip. Typically, the early dip has not been observed consistently in BOLD fMRI but is a robust feature in OIS. The source of this discrepancy is also unknown. Additional systematic studies are necessary to understand the differences between BOLD and OIS signal sources at a submillimeter scale.
Footnotes
-
This work was supported by National Institutes of Health (NIH) Grants NS44589, EB003324, and EB003375 and by the McKnight Foundation. The 9.4 T MR system was funded in part by NIH Grant RR17239. We thank Ping Wang and Michelle Tasker for animal preparation and Kristy Hendrich for MR system management.
- Correspondence should be addressed to Seong-Gi Kim, Department of Radiology, University of Pittsburgh, 3025 East Carson Street, Pittsburgh, PA 15203. kimsg{at}pitt.edu