## Abstract

Neuroimaging studies of functional magnetic resonance imaging (fMRI) and electrophysiology provide the linkage between neural activity and the blood oxygenation level–dependent (BOLD) response. Here, BOLD responses to light flashes were imaged at 11.7T and compared with neural recordings from superior colliculus (SC) and primary visual cortex (V1) in rat brain—regions with different basal blood flow and energy demand. Our goal was to assess neurovascular coupling in V1 and SC as reflected by temporal/spatial variances of impulse response functions (IRFs) and assess, if any, implications for general linear modeling (GLM) of BOLD responses. Light flashes induced high magnitude neural/BOLD responses reproducibly from both regions. However, neural/BOLD responses from SC and V1 were markedly different. SC signals followed the boxcar shape of the stimulation paradigm at all flash rates, whereas V1 signals were characterized by onset/offset transients that exhibited different flash rate dependencies. We find that IRFSC is generally time-invariant across wider flash rate range compared with IRFV1, whereas IRFSC and IRFV1 are both space invariant. These results illustrate the importance of measured neural signals for interpretation of fMRI by showing that GLM of BOLD responses may lead to misinterpretation of neural activity in some cases.

## Introduction

Functional magnetic resonance imaging (fMRI) is widely used to map brain activity in humans and animals because of its high spatiotemporal resolution and extensive volumetric coverage. The popular use of fMRI relies on the fact that the blood oxygenation level–dependent (BOLD) signal is a noninvasive, albeit indirect, indicator of neural activity (Ogawa et al. 1990). The BOLD signal is actually an index of the hyperemic response, which has both hemodynamic (blood flow and volume) and metabolic (oxidative energy metabolism) contributions (Hyder et al. 2001). Hence the interpretation of fMRI results depends on the connection between the hyperemic response and neural activity (for a recent review, see Hyder et al. 2010), as indicated by neural spiking rate of multiunit activity (MUA) and/or integrative processes reflected by local field potentials (LFPs).

In experimental settings where neural activity cannot be readily measured in conjunction with fMRI, general linear modeling (GLM) of the BOLD response is used to decipher the underlying, albeit unmeasured, neural activity (Friston 1998). Today, GLM is the most prevalent method to examine fMRI time series from human experiments, both with and without task (Frackowiak 2004; Huettel et al. 2004). When applied in the most common context of activation detection, the GLM consists of an assumed underlying “neural activity function” (typically where the stimulus paradigm itself is used as the neural time course) that is convolved with a hemodynamic response function (HRF) representing neurovascular coupling. Coupling linearity is explicitly assumed, and hence, in its most commonly applied form, GLM is undefined for nonlinear effects, whether they are temporal or spatial in nature.

Studies exploring the neurophysiological basis of fMRI—now being conducted in many laboratories for various systems and species (Logothetis et al. 2001; Smith et al. 2002; Kim et al. 2004; Maandag et al. 2007; van der Linden et al. 2007; Goense and Logothetis 2008; Huttunen et al. 2008; Pawela et al. 2008; Liu et al. 2009)—aim to describe the link between neural activity and the BOLD response. When both electrical and vascular signals are measured in identical stimulation conditions, one may quantify this link in the form of an impulse response function (IRF). An IRF, generated by deconvolution of the measured neural activity from the measured BOLD response, characterizes the neurovascular coupling for a given region (Iadecola 2004) and can even be used to calibrate the BOLD response (Herman et al. 2009; Sanganahalli et al. 2009b). Another important IRF research direction for fMRI, however, is to determine how widely applicable is the suggested link between changes in neural activity and the BOLD response across experimental conditions (e.g., across regions or variable stimuli), because IRF estimation is an identification method applicable to linear systems (Hespanha 2009). To accurately interpret fMRI data, it is vital to know whether the time-invariance assumption holds for a wide range of stimulus parameters and how space-invariant IRFs are across regions with divergent basal hemodynamics and metabolic demand (Boynton et al. 1996; Engel et al. 1997).

Based on the hypotheses that neurovascular couplings in the rat visual system are similar 1) across brain regions and 2) for varying flash rates, we sought to evaluate temporal and spatial variances of IRFs and their implications for GLM of BOLD responses. We measured BOLD responses to light flashes at 11.7T in rat brain and compared these signals with extracellular recordings from superior colliculus (SC) and primary visual cortex (V1), two regions which have different basal blood flow and energy demand (Iadecola et al. 1983; Klein et al. 1986; Borowsky and Collins 1989; Maeda et al. 1998). A binocular stimulus delivery system allowed variations of light color and flash rate to evoke reproducible BOLD and neural responses at a single-trial level. Although neural and BOLD responses to blue/green light generally agreed with each other, responses from V1 and SC were noticeably different. We find that time-invariance holds for a wide range of flash rates in SC but not in V1, whereas IRFSC and IRFV1 are both space invariant. These results suggest caution in using GLM of BOLD responses because of the potential to misconstrue underlying neural activity for varying stimulus rates.

## Materials and Methods

### Animal Preparation

All procedures were performed in accordance with protocols approved by the Yale University Institutional Animal Care and Use Committee, in agreement with the National Institutes of Health Guide for the Care and Use of Laboratory Animals. Experiments were conducted on artificially ventilated (70% N2O/30% O2) adult male rats (Long Evans; 220–350 g; Charles River, Wilmington, MA). A femoral artery was cannulated (PE-50) for monitoring physiological parameters (pCO2, pO2, pH, and blood pressure), and a femoral vein was cannulated (PE-10) to administer d-tubocurarine chloride (initial 0.5 mg/kg; supplemental 0.25 mg/kg/h) and other substances. The rats were initially anesthetized with isoflurane (3–4%) to minimize stress caused by injection of urethane (1.3 g/kg, intraperitoneal), which was used to maintain anesthesia throughout the experiment. Additional doses of urethane (0.13 g/kg, intravenous) were administered if necessary depending on the duration of the experiment. Ventilation parameters were adjusted to maintain normal physiology. The animal core temperature was measured by a rectal probe and controlled by a heated water blanket maintained at ∼37 °C. Adequate analgesia was ensured by monitoring the pain response to an automated electrical (5 mA, 0.3 ms, 3 Hz, 1–5 s) tail pinch every 15 min in general, but not during the visual stimulation protocol.

Details of the stimulus delivery system are given elsewhere (Sanganahalli, Bailey, et al. 2009), but briefly summarized here and illustrated in Supplementary Figure 1. A computer-controlled CED μ1401 analog-to-digital converter unit (CED unit) and Spike2 software (Cambridge Electronic Design, Cambridge, UK) were used to control stimulator output and record measured data for off-line analysis. The design of the setup allows control of light intensity and duration of flashes continuously, whereas wavelength (color) could be varied between trials by selecting a different light source amongst the chosen light emitting diodes (LEDs; Luxeon Star III, Lumileds Lighting, LLC, San Jose, CA; wavelengths blue = 455 nm, green = 530 nm, and red = 627 nm). We used acrylic collimator lenses to target the wide-angle LED light into optical cables (Fiber Optic Products, Clearlake Oaks, CA), which attached to cradles in front of the LEDs. Two thin (inner/outer diameters of 1.0/2.3 mm) polyethylene-coated fibers for independent stimulation of each eye were led into the scanner room (Faraday cage) through small access holes in the sidewall (back wall). A fiberglass platform in front of the animal held in place and allowed the positioning of the thin cables by means of cylindrical anchoring constructions, which also slightly dispersed the light thus increasing the viewing angle of the stimuli. We thus estimated the portion of the visual field subtended by each stimulus to be approximately 35–50° azimuth and 10–40° elevation relative to the horizontal meridian. However, as each emitted beam of light illuminated the entire eyeball, the stimuli essentially cover the full visual field. The illuminance of each color stimulus was measured before each experiment, in order to ensure equal intensities for left and right visual fields. The calibration was performed in Lux units using a standard light meter (Extech Instruments, Waltham, MA). We emphasize that we did not strive to attain equal illuminances between colors, as the Lux is a measure dependent on human subjective intensity ratings, and thereby ultimately on the properties of the human retina. The dim ambient lighting conditions, that is, only reflections from distant halogen ceiling lights, in the MR as well as electrophysiology settings, were not considered to bias the detection of the stimuli toward any particular wavelength. Beyond these low-light settings, the rats were not dark adapted for the measurements. A standard block stimulation protocol was used in all experiments with at least 300 s between consecutive trials.

### Electrophysiology Experiments

Rats were mounted on a stereotaxic frame (David Kopf Instruments, Tujunga, CA) placed on a vibration-free table inside a Faraday cage. The scalp and the “galea aponeurotica” were removed and small burr holes were drilled for insertion of high-impedance tungsten microelectrodes (2–4 MΩ; FHC, Bowdoinham, ME) that measured neural electrical signals. The coordinates for electrode placement were guided by our fMRI data and the atlas of Paxinos and Watson (1998). With reference to bregma and the sagittal midline plane, the electrodes were placed in the following coordinates: V1 (−6.5 mm posterior, 3.5 mm lateral, ∼1 mm from dura) and SC (−7.0 mm posterior, 0.8 mm lateral, 2.4–2.6 mm from dura). Electrical signals were digitized at 20 kHz with CED μ1401 using Spike2. LFP and MUA were obtained by splitting the electrical signal into low (<300 Hz) and high (400 Hz–10 kHz) frequency bands using an online analogue Butterworth filter (24 dB/oct attenuation).

Off-line data analysis was performed using custom scripts in Matlab (The MathWorks, Inc., USA). The MUA signal was squared and resampled to a sampling rate of 8 Hz using the Matlab function decimate. The resampling procedure was recursive in several steps, each decimation stage being at maximum 10. The function decimate was applied using a 30th order FIR lowpass filter and forward and backward filtering passes to avoid distorting the phase of the signals (i.e, because the timing of peaks would otherwise be affected). Then finally, the signals were smoothed using a 1.25-s zero-phase moving average FIR filter, and square rooted to obtain a root-mean-square (RMS) representation of the signal power similar to Stark and Abeles (2007).

### fMRI Experiments

The fMRI data were obtained on a modified 11.7T Bruker horizontal-bore spectrometer (Bruker AVANCE, Billerica, MA) using a 1H surface coil (diameter 1.4 cm). Details of fMRI measurements are discussed elsewhere (Sanganahalli et al. 2008). Briefly, BOLD contrast images were acquired using echo-planar imaging with sequential sampling (Hyder et al. 1995) and the following parameters: repetition time = 1000 ms; echo time = 15 ms; field of view = 2.56 × 2.56 cm2; image matrix = 64 × 64; number of slices = 5; slice thickness = 1 mm; and flip angle = 60°. The neuroanatomy was imaged with the rapid acquisition relaxation enhanced (Hennig et al. 1986) or the fast low-angle single shot (Frahm et al. 1986) contrast sequences. In the anesthetized rat, different gradient amplitude and pulsing rates do not change the blood pressure and electrical recordings outside the magnet (on the bench) and inside the magnet (during an fMRI scan) are similar (Englot et al. 2008).

Analysis of fMRI data was performed on a trial-by-trial basis using custom-written routines in Matlab. Trials with excessive motion artifacts, defined by translational movement in excess of a quarter of the voxel size and estimated using center-of-mass analysis, were rejected. Each voxel's BOLD data were subjected to a Student's t-test comparing the mean of the 30 s baseline to that of the stimulation period. A significance threshold of P < 0.01 (uncorrected for multiple comparisons) was used to create maps showing the activated tissue area. Thresholded t-maps were overlaid on corresponding structural images for visualization and used to define a region-of-interest (ROI) for extraction of the BOLD signal time course. For responses lacking sustained stimulus-induced signal increases, signal time courses were extracted using an ROI defined by other stimuli. Before averaging across activated voxels, the data were normalized to the baseline period, thus yielding plots of ΔS/S.

### Estimation of IRF

We modeled the 1-Hz BOLD signal as the output with MUA signal as the input.

(1)
$BOLDiR(t)≡hR(t)⊗MUAjR(t),$
where $BOLDiR(MUAjR)$ is the ith (jth) BOLD (MUA) records in region R:{V1, SC}, the transfer function hR(t) is independent of i and j, and ⊗ denotes convolution. Since the general deconvolution problem of estimating h(t) from BOLD-MUA pairs is ill-conditioned, we used the gamma variate function (Madsen 1992) as the IRF
(2)
$h(t)=A(ttpe1ttp)α,$
where A is a scaling constant, tp is the peak time, and α is related to the width of the function. We used standard nonlinear optimization (Matlab function lsqcurvefit with trust region reflective algorithm) to find the parameters {A,tp,α} that minimize the squared difference between measured $(BOLDiR(t))$ and modeled $BOLD(hR(t)⊗MUAjR(t))$ for i,j = {1,. . . ,11} (121 pairs). The final parameter estimates were obtained by averaging over all fits after regularizing against outliers by removing the 5% of the largest and smallest tp × α-products (corresponding to late sharp and early broad peaks, respectively).

### GLM of BOLD Response

A model-based analysis was performed on a subset of the data to examine frequency-dependent effects. A single gamma variate was chosen as the canonical HRF with tp =3.23 s and α = 11.33. The boxcar regressors of the GLM were defined as follows and subsequently convolved with the canonical HRF: initial onset response (0–2 s), sustained response (0–30 s), offset response (31–33 s), and poststimulus signal elevation (31–90 s). Parameter estimates for ROI time courses were obtained using ordinary least squares in Matlab. Statistical parametric maps were calculated voxel-by-voxel using NeuroLens (www.neurolens.org) after smoothing with a 0.8 mm full-width at half-maximum Gaussian spatial filter.

## Results

### Color Sensitive Responses at a Single-Trial Level

Strong BOLD responses were observed in V1 and SC. Figure 1A illustrates the activation maps above the selected threshold (30 s stimulation vs. 30 s baseline; P < 0.01) for 2 consecutive trials of 1-Hz visual flash stimulation in one animal. Cortical responses extended into surrounding secondary visual processing areas (V2), but peaked in slices approximately 6–7 mm posterior to bregma, ∼3.5 mm lateral to midline, and at a depth of 1–1.2 mm relative to dura, which corresponds approximately to layers IV–V in the V1 area. Retinal ganglion cells project directly to the most dorsal layers of the SC, where maximal BOLD responses were observed, approximately 7–8 mm posterior to bregma, 0–1 mm lateral to midline, and at a depth of 2.2–2.8 mm relative to dura. These coordinates were used for subsequent microelectrode examinations of neural responses to identical stimulation conditions. We limited our focus on V1 and SC responses because the surface transceiver coil was positioned optimally to detect BOLD signals from these regions. Thus low signal-to-noise ratio was observed in the more anterior dorsal lateral geniculate (DLG) nucleus of the thalamus (∼4–5 mm posterior to bregma).

Figure 1.

BOLD and neural responses to light stimulation in the rat visual system. (A) BOLD activation maps from 2 trials of green light 1 Hz flash stimulation executed on the same animal (two-sample Student's t-test between 30 s baseline and 30 s stimulation periods, thresholded at P = 0.01, i.e., −log(P) = 4.6). Significant signal increases were observed in the middle cortical layers of the primary (V1) and secondary (V2) visual areas, and in the dorsal layers of the superior colliculus (SC). The z-coordinates between slices are relative to bregma (Paxinos and Watson 1998). (B) Single trial time courses of BOLD (top) and neural (bottom) responses to red, green, and blue light (1-Hz stimulation blocks). BOLD and MUA data from V1 (right) and SC (left) show color dependence, where close resemblance of BOLD and neural responses are observed with green and blue light. In comparison, red light induces weaker responses in both V1 and SC. The black horizontal bar below each plot marks the stimulation period (30 s). For illustration of the custom-made nondark adapted binocular stimulus delivery system and regional overlap of BOLD activation maps with different light colors, see Supplementary Figures 1 and 2, respectively.

Figure 1.

BOLD and neural responses to light stimulation in the rat visual system. (A) BOLD activation maps from 2 trials of green light 1 Hz flash stimulation executed on the same animal (two-sample Student's t-test between 30 s baseline and 30 s stimulation periods, thresholded at P = 0.01, i.e., −log(P) = 4.6). Significant signal increases were observed in the middle cortical layers of the primary (V1) and secondary (V2) visual areas, and in the dorsal layers of the superior colliculus (SC). The z-coordinates between slices are relative to bregma (Paxinos and Watson 1998). (B) Single trial time courses of BOLD (top) and neural (bottom) responses to red, green, and blue light (1-Hz stimulation blocks). BOLD and MUA data from V1 (right) and SC (left) show color dependence, where close resemblance of BOLD and neural responses are observed with green and blue light. In comparison, red light induces weaker responses in both V1 and SC. The black horizontal bar below each plot marks the stimulation period (30 s). For illustration of the custom-made nondark adapted binocular stimulus delivery system and regional overlap of BOLD activation maps with different light colors, see Supplementary Figures 1 and 2, respectively.

We anticipated, based on known retinal morphology of the rat (see Discussion), that responses to red light flashes would be poorer in eliciting neural and BOLD responses compared with green and/or blue light flashes. Indeed, Figure 1B illustrates that our paradigm is sufficiently sensitive to detect color-dependent neural and BOLD responses at a single-trial level. Figure 1B shows the effect of varying stimulus color (wavelength) for 2 individual animals (one for BOLD, another for MUA). Blue and green light consistently produced responses of similar magnitude and temporal evolution in both V1 and SC. Cortical responses to red light tend to be slightly weaker and less robust compared with other colors, particularly in the SC (Figure 1B, left). The spatial overlap of the BOLD responses to different colors is practically complete, with fewer significantly modulated voxels during red than blue/green stimulation, particularly in the SC (Supplementary Fig. 2). We, therefore, analyzed trials with green and blue light, whereas red light trials are not included in subsequent analysis.

### Differential V1 and SC Response Characteristics

Because V1 and SC signals were measured simultaneously, whether they are in fMRI studies with the rat inside the magnet or electrophysiology recordings with the rat on the bench, we quantified BOLD and neural response reproducibility over trials and animals by calculating the pairwise cross correlations of V1 and SC time courses in 11 representative trials in 6 animals and plotted the squared correlation coefficient (Pearson R2) value for each V1 × V1, SC × SC, and V1 × SC trial pair (gray scale coded 11 × 11 matrix). We used the 1-Hz flash rate data for this analysis. The darker shade in the correlation matrices in Figure 2A indicate that stronger correlations are observed within than across regions, more so for BOLD data than MUA data (presumably because the initial onset response of the former occupies a larger portion of the 30-s stimulation period) used in the cross correlations. The mean (±1 standard deviation [SD]) R2 between individual trials is 0.65 (±0.11) for V1 BOLD, 0.65 (±0.12) for SC BOLD, 0.76 (±0.09) for V1 MUA and 0.87 (±0.07) for SC MUA. In comparison, the cross regional (V1 × SC) trial-by-trial correlations are smaller for BOLD (0.37 ± 0.15) and MUA (0.68 ± 0.08) data (see bar graph insets in Fig. 2A).

Figure 2.

BOLD and neural temporal response profiles in V1 and SC differ significantly. (A) Trial-by-trial reproducibility analysis of 11 examples of BOLD and neural responses to 1-Hz flash stimulation (data from 6 animals). Each matrix, BOLD signal left, and MUA right, illustrates the R2-values for pairwise cross correlations between single trial time courses from the same region (V1 × V1 and SC × SC, i.e., intraregional) as well as those between time courses from different regions (V1 × SC, i.e., interregional). The adjacent bar graphs summarize these results, indicating higher correlations within a given region than between regions, due to different signal time courses in V1 and SC. (B and C) Average time courses (±SD) of BOLD signal and MUA measured from V1 (black) and SC (gray) for 1-Hz flash stimulation. The time courses of BOLD and MUA from V1 show a prominent onset peak, and a weaker but noticeable offset transient (arrows). The gray horizontal bar represents the stimulation period (30 s). For illustration of the corresponding LFP data, see Supplementary Figure 3.

Figure 2.

BOLD and neural temporal response profiles in V1 and SC differ significantly. (A) Trial-by-trial reproducibility analysis of 11 examples of BOLD and neural responses to 1-Hz flash stimulation (data from 6 animals). Each matrix, BOLD signal left, and MUA right, illustrates the R2-values for pairwise cross correlations between single trial time courses from the same region (V1 × V1 and SC × SC, i.e., intraregional) as well as those between time courses from different regions (V1 × SC, i.e., interregional). The adjacent bar graphs summarize these results, indicating higher correlations within a given region than between regions, due to different signal time courses in V1 and SC. (B and C) Average time courses (±SD) of BOLD signal and MUA measured from V1 (black) and SC (gray) for 1-Hz flash stimulation. The time courses of BOLD and MUA from V1 show a prominent onset peak, and a weaker but noticeable offset transient (arrows). The gray horizontal bar represents the stimulation period (30 s). For illustration of the corresponding LFP data, see Supplementary Figure 3.

The average V1 and SC response time courses to the 1-Hz flash rate are shown for both BOLD (Fig. 2B) and MUA (Fig. 2C). The V1 neural recordings exhibit a sharp onset transient (tpeak = 327 ± 79 ms), an initial rapid decay (peak-to-plateau = 791 ± 158 ms), followed by a plateau of more gradual decline in response amplitude until cessation of the stimulus. The V1 mean plateau-to-peak ratio is 38.0 ± 12.0%. In contrast, the SC shows a rapid rise to peak (382 ± 275 ms) but only moderate decline during sustained stimulation. Similar dynamics were obtained for the 80–200 Hz band of the LFP signal (for a full time-frequency decomposition, see Supplementary Fig. 3 and for details on the wavelet filtering applied, see Supplementary Methods). The BOLD data are in excellent agreement with the neural signals, indicating similar neurovascular coupling mechanisms in the 2 regions during 1 Hz stimulation. The V1 and SC BOLD signal peaks occur at 5.2 ± 0.7 s and 5.8 ± 1.4 s, respectively. The V1 plateau-to-peak ratio is 34.0 ± 11.0%. Note that the V1 MUA signal, and correspondingly the V1 BOLD, exhibits a transient event to stimulus train offset (30–32 s for MUA; see arrows in Fig. 2B,C).

The averaged BOLD and neural temporal response profiles in response to higher flash rates (5 and 10 Hz, respectively) are shown in Figure 3A,B. Several features of the response dynamics remain similar to the 1 Hz flash rate data (Fig. 2), but differences are also evident at both 5- and 10-Hz flash rates for SC and V1. A remarkable feature is that the plateau of the BOLD response in V1 (but not its neural response) is strongly attenuated at the 10 Hz flash rate, revealing neurovascular uncoupling (i.e., MUA > BOLD during sustained response in V1). Likewise in SC, the prominent offset BOLD response (but not in its neural response) is augmented at the 10-Hz flash rate, revealing neurovascular uncoupling but in the opposite direction from that observed in V1 (i.e., BOLD > MUA during offset response in SC). The regional neurovascular uncouplings for onset, sustained, and offset parts of the responses at 10-Hz flash rate are summarized in Table 1. For both of these features observed in SC and V1, it can be argued that the neurovascular uncoupling is maximum at 10-Hz flash rate, but the trend is partly detectable at 5 Hz (but absent at 1 Hz) flash rate. Supplementary Figures 4 and 5 represent the time-frequency decompositions of the 5- and 10-Hz LFP data and illustrate the good agreement between the MUA data shown in Figure 3 and the 80–200 Hz band LFP signals.

Table 1

Summary of neurovascular coupling and uncoupling observed during 10-Hz flash rate

 V1 SC Onset (∼3 s in duration) Coupling observed Coupling observed Sustained (>25 s in duration) Coupling NOT observed Coupling observed Offset (∼3 s in duration) Coupling observed Coupling NOT observed
 V1 SC Onset (∼3 s in duration) Coupling observed Coupling observed Sustained (>25 s in duration) Coupling NOT observed Coupling observed Offset (∼3 s in duration) Coupling observed Coupling NOT observed
Figure 3.

Higher flash rates induce variable BOLD and neural temporal response profiles in V1 and SC. Average time courses (±SD) of BOLD signal (top) and MUA (bottom) measured from V1 (black) and SC (gray) for (A) 5- and (B) 10-Hz flash rate stimulation (data from 6 animals). Data obtained at the higher flash rates of 5 Hz and 10 Hz demonstrate different dynamics from data obtained at 1 Hz. First, there are prominent offset responses from V1 in the BOLD and MUA data. Second, the plateau BOLD response from V1 is significantly attenuated at the10-Hz flash rate. Third, the BOLD and MUA data from SC do not show significant onset responses, but there are some offset responses at higher flash rates. For summary of neurovascular couplings and uncouplings for different flash rates in V1 and SC, see Table 1. The gray horizontal bar represents the stimulation period (30 s). For illustration of the corresponding LFP data, see Supplementary Figures 4 and 5.

Figure 3.

Higher flash rates induce variable BOLD and neural temporal response profiles in V1 and SC. Average time courses (±SD) of BOLD signal (top) and MUA (bottom) measured from V1 (black) and SC (gray) for (A) 5- and (B) 10-Hz flash rate stimulation (data from 6 animals). Data obtained at the higher flash rates of 5 Hz and 10 Hz demonstrate different dynamics from data obtained at 1 Hz. First, there are prominent offset responses from V1 in the BOLD and MUA data. Second, the plateau BOLD response from V1 is significantly attenuated at the10-Hz flash rate. Third, the BOLD and MUA data from SC do not show significant onset responses, but there are some offset responses at higher flash rates. For summary of neurovascular couplings and uncouplings for different flash rates in V1 and SC, see Table 1. The gray horizontal bar represents the stimulation period (30 s). For illustration of the corresponding LFP data, see Supplementary Figures 4 and 5.

### IRFV1 and IRFSC for 1-Hz Flash Rate

To quantify the extent of neurovascular coupling in V1 and SC, we estimated the parameters of the IRF (i.e., h(t) in equations (1) and (2)) that lead to optimal fits between the measured and modeled BOLD responses based on the convolution between MUA and the IRF. We used the 1 Hz data shown in Figure 2 for the estimation of h(t), leading to 121 pairs of MUA and BOLD signals in both V1 and SC. For numerical tractability, we limited the optimization to the 3 parameters (amplitude, peak, and width) of the gamma variate family of functions. Furthermore, the fitting time of each data set was limited to [−5, 40] s (i.e., to mainly the stimulation period; for details on the estimation procedure, see Materials and Methods).

The estimated IRFs for V1 and SC data at 1-Hz flash rate are shown in Figure 4. The peak and width parameters of the IRFs are of similar magnitude and estimates statistically indistinguishable (overlapping SDs): IRFV1 and IRFSC peak at 2.5 ± 0.5 s and 2.0 ± 1.1 s, respectively, with widths at half height of 1.5 ± 0.4 s and 1.4 ± 0.3 s. These large variations of the peaks and widths of the IRFs are due to intersubject variations. Thus, we conclude that there is no evidence supporting different neurovascular (BOLD-MUA) couplings in V1 and SC during 1-Hz flash stimulation. We repeated the estimation procedure using the power in the 80–200 Hz band of the LFP signal as the neural input function (for equations (1) and (2)). The IRFs generated for the “high-gamma” LFP and BOLD responses are in close agreement to those obtained for BOLD and MUA responses (data not shown). We also calculated IRFs for SC and V1 at 5-Hz flash rate and found negligible differences from IRFs with 1-Hz flash rate (data not shown). However (linear) IRFs at 10-Hz flash rate are not defined due to the inherent neurovascular uncouplings observed in V1 during the sustained response (MUA > BOLD) and in SC during the offset response (BOLD > MUA) (see Table 1 and Fig. 3).

Figure 4.

Average IRFs for V1 and SC (i.e., IRFV1 and IRFSC represented by solid and dashed lines, respectively) obtained by deconvolution between single-trial BOLD and MUA data measured during 1 Hz stimulation (data from 6 animals). See text for details on the gamma variate function. IRFV1 and IRFSC peak at 2.5 ± 0.5 s and 2.0 ± 1.1 s, respectively, with widths at half height of 1.5 ± 0.4 s and 1.4 ± 0.3 s. The gray shading represents the SD of the estimated IRFs.

Figure 4.

Average IRFs for V1 and SC (i.e., IRFV1 and IRFSC represented by solid and dashed lines, respectively) obtained by deconvolution between single-trial BOLD and MUA data measured during 1 Hz stimulation (data from 6 animals). See text for details on the gamma variate function. IRFV1 and IRFSC peak at 2.5 ± 0.5 s and 2.0 ± 1.1 s, respectively, with widths at half height of 1.5 ± 0.4 s and 1.4 ± 0.3 s. The gray shading represents the SD of the estimated IRFs.

### Temporal Invariance of the IRF

We tested the temporal invariance of the IRFs estimated for the 1-Hz flash rate (Fig. 4) by applying them to data collected during 5-Hz and 10-Hz flash rates. Thus, we convolved the MUA signals recorded during stimulation at these higher flash rates with the region-specific IRFs shown in Figure 4 (i.e., IRFV1 for V1 data and IRFSC for SC data) to provide comparisons of modeled and measured BOLD responses (Fig. 5). Assumption of linearity (and temporal invariance) holds if the modeled BOLD responses can be scaled to match the measured BOLD responses at all flash rates. Comparisons of MUA and BOLD responses (modeled and measured single trials) are shown for 1-, 5-, and 10-Hz flash rates in Figures 5A,B,C, respectively. The predictions for SC are accurate within experimental noise for all flash rates, whereas those for V1 concur with measurements only for 1- and 5-Hz flash rates, but not for flash rate of 10 Hz. Therefore, the IRFSC (in Fig. 4) is time invariant at all applied flash rates, whereas the IRFV1 (Fig. 4) is time invariant except at the 10-Hz flash rate. We repeated the analysis shown in Figure 5 for single trials on the whole data set, and plotted the modeled versus measured BOLD responses in the period from −5 s to +40 s. Consistent with the above observations, the data points for both SC and V1 comparisons lie close to the line of identity, except for the case of 10-Hz flash rate stimulus in V1, where the prediction consistently overestimates the measured BOLD signal (Supplementary Fig. 6A,B).

Figure 5.

Comparison between measured and modeled BOLD responses to test the temporal invariance of modeled SC and V1 signals (i.e., IRFSC and IRFV1). Data for V1 and SC are shown in the right and left columns, respectively. Time courses of measured BOLD responses (lower panel, gray shade) for (A) 1-, (B) 5-, and (C) 10-Hz flash rates are juxtaposed against the time courses of modeled BOLD responses (lower panel, black trace), which were generated by convolution between the respective MUA responses (upper panel, black trace) and corresponding IRFs for 1 Hz stimulation for that region (Fig. 4). In other words, the modeled BOLD responses from V1 were obtained using IRFV1 and the modeled BOLD responses from SC were obtained using IRFSC. Good agreement between modeled and measured BOLD responses is found in SC at all flash rates examined, whereas only the 1 and 5 Hz responses are adequately modeled by IRFV1. See Supplementary Figure 6A,B for illustration of correspondence of measured and modeled signals from all animals. The black horizontal bar represents the stimulation period (30 s). Data from a single representative animal.

Figure 5.

Comparison between measured and modeled BOLD responses to test the temporal invariance of modeled SC and V1 signals (i.e., IRFSC and IRFV1). Data for V1 and SC are shown in the right and left columns, respectively. Time courses of measured BOLD responses (lower panel, gray shade) for (A) 1-, (B) 5-, and (C) 10-Hz flash rates are juxtaposed against the time courses of modeled BOLD responses (lower panel, black trace), which were generated by convolution between the respective MUA responses (upper panel, black trace) and corresponding IRFs for 1 Hz stimulation for that region (Fig. 4). In other words, the modeled BOLD responses from V1 were obtained using IRFV1 and the modeled BOLD responses from SC were obtained using IRFSC. Good agreement between modeled and measured BOLD responses is found in SC at all flash rates examined, whereas only the 1 and 5 Hz responses are adequately modeled by IRFV1. See Supplementary Figure 6A,B for illustration of correspondence of measured and modeled signals from all animals. The black horizontal bar represents the stimulation period (30 s). Data from a single representative animal.

### Spatial Invariance of the IRF

The demonstration of space-invariant neurovascular coupling is shown in Figure 6, in which IRFV1 is applied on the data obtained form the SC and IRFSC on the data obtained from V1. Assumption of linearity (and spatial invariance) holds if the modeled BOLD responses can be scaled to match the measured BOLD responses from both regions. As can be seen in Figure 6A,B,C, respectively, for 1-, 5-, and 10-Hz flash rates, the effect of interchanging the IRF of one region with that of the other is minimal and within experimental error at all flash rates (compare Fig. 5 with Fig. 6). We repeated the analysis on all trials and plotted modeled BOLD responses against measured BOLD responses (Supplementary Fig. 6C,D). We conclude that IRFSC and IRFV1 (Fig. 4) are space invariant within experimental error.

Figure 6.

Comparison between measured and modeled BOLD responses to test the spatial invariance of IRFSC and IRFV1. Data for V1 and SC are shown in the right and left columns, respectively. Time courses of measured BOLD responses (lower panel, gray shade) for (A) 1-, (B) 5-, and (C) 10-Hz flash rates are juxtaposed against the time courses of modeled BOLD responses (lower panel, black trace), which were generated by convolution between the respective MUA responses (upper panel, black trace) and corresponding IRFs for 1 Hz stimulation for that region (Fig. 4). In other words, the modeled BOLD responses from V1 were obtained using IRFSC and the modeled BOLD responses from SC were obtained using IRFV1. See Supplementary Figure 6C,D for illustration of space-invariant IRF analysis from all animals. The black horizontal bar represents the stimulation period (30 s). Data from a different representative animal from that in Figure 5.

Figure 6.

Comparison between measured and modeled BOLD responses to test the spatial invariance of IRFSC and IRFV1. Data for V1 and SC are shown in the right and left columns, respectively. Time courses of measured BOLD responses (lower panel, gray shade) for (A) 1-, (B) 5-, and (C) 10-Hz flash rates are juxtaposed against the time courses of modeled BOLD responses (lower panel, black trace), which were generated by convolution between the respective MUA responses (upper panel, black trace) and corresponding IRFs for 1 Hz stimulation for that region (Fig. 4). In other words, the modeled BOLD responses from V1 were obtained using IRFSC and the modeled BOLD responses from SC were obtained using IRFV1. See Supplementary Figure 6C,D for illustration of space-invariant IRF analysis from all animals. The black horizontal bar represents the stimulation period (30 s). Data from a different representative animal from that in Figure 5.

### GLM of the BOLD Response

We constructed a GLM for the rat visual system that embodies the main features observed in V1 and SC BOLD responses: onset (0–2 s) and offset (31–33 s) transients, a sustained component during stimulation (0–30 s), and a poststimulus baseline shift (31–90 s). Each of these 4 components of the model was convolved with a gamma variate-based canonical HRF using parameters set to the average of IRFV1 and IRFSC (tp = 3.23 s, α = 11.33). Results for a single representative animal are shown in Figure 7. An identical analysis was performed for 5 other data sets, and the group results are summarized in Table 2.

Table 2

Values of parameter estimates of the four phases GLM following 1-, 5-, and 10-Hz flash rate stimulation (n = 6)

 V1 (% ± SD) SC (% ± SD) 1 Hz 5 Hz 10 Hz 1 Hz 5 Hz 10 Hz Onset 2.01 (±0.59) 0.95 (±0.33) 0.76 (±0.46) 0.59 (±0.70) 0.46 (±0.64) −0.18 (±0.42) Sustained 0.79 (±0.19) 0.43 (±0.24) −0.05 (±0.47) 1.66 (±0.29) 1.80 (±0.40) 1.79 (±0.45) Offset 0.69 (±0.29) 1.78 (±0.44) 2.25 (±0.78) 0.51 (±0.60) 0.77 (±0.52) 0.99 (±0.66) Poststimulation −0.16 (±0.12) 0.01 (±0.20) −0.19 (±0.23) −0.07 (±0.29) 0.70 (±0.50) 1.06 (±0.36)
 V1 (% ± SD) SC (% ± SD) 1 Hz 5 Hz 10 Hz 1 Hz 5 Hz 10 Hz Onset 2.01 (±0.59) 0.95 (±0.33) 0.76 (±0.46) 0.59 (±0.70) 0.46 (±0.64) −0.18 (±0.42) Sustained 0.79 (±0.19) 0.43 (±0.24) −0.05 (±0.47) 1.66 (±0.29) 1.80 (±0.40) 1.79 (±0.45) Offset 0.69 (±0.29) 1.78 (±0.44) 2.25 (±0.78) 0.51 (±0.60) 0.77 (±0.52) 0.99 (±0.66) Poststimulation −0.16 (±0.12) 0.01 (±0.20) −0.19 (±0.23) −0.07 (±0.29) 0.70 (±0.50) 1.06 (±0.36)
Figure 7.

Comparison between experimentally measured BOLD responses and modeled BOLD responses created from GLM. Time courses (right column) of measured BOLD responses (pink shade) from V1 (top trace) and SC (bottom trace) for (A) 1-, (B) 5-, and (C) 10-Hz flash rates are juxtaposed with the BOLD activation maps (left column) for the 4 phases of the response: onset, sustained, offset, and poststimulation. The modeled BOLD responses were generated by summing the time courses generated from GLM of the 4 components, where the canonical HRF used for GLM was the average of IRFSC and IRFV1 (in Fig. 4), but very similar results are obtained if either of IRFSC and IRFV1 is used (data not shown). The gray horizontal bar represents the stimulation period (30 s). Data from a single representative animal; for summary statistics for 6 animals, see Table 2. poststimulation.

Figure 7.

Comparison between experimentally measured BOLD responses and modeled BOLD responses created from GLM. Time courses (right column) of measured BOLD responses (pink shade) from V1 (top trace) and SC (bottom trace) for (A) 1-, (B) 5-, and (C) 10-Hz flash rates are juxtaposed with the BOLD activation maps (left column) for the 4 phases of the response: onset, sustained, offset, and poststimulation. The modeled BOLD responses were generated by summing the time courses generated from GLM of the 4 components, where the canonical HRF used for GLM was the average of IRFSC and IRFV1 (in Fig. 4), but very similar results are obtained if either of IRFSC and IRFV1 is used (data not shown). The gray horizontal bar represents the stimulation period (30 s). Data from a single representative animal; for summary statistics for 6 animals, see Table 2. poststimulation.

The left column of Figure 7 illustrates thresholded statistical parametric maps (left panel) for each of the 4 model components for 1- (Fig.7A), 5- (Fig. 7B), and 10-Hz (Fig. 7C) flash rates, respectively. The first observation is that the V1 transient and the 2 sustained components (when present) are colocalized, that is, all 3 response components are generated by the same region of V1 and presumably in the same neuronal population from which the electrical signals (i.e., MUA and LFP) were recorded. Compared with the sustained SC response, the poststimulus signal level elevation observed in SC is more restricted to midline structures. The right column of Figure 7 demonstrates the GLM fits for each of the 4 model components for 1- (Fig. 7A), 5- (Fig. 7B), and 10-Hz (Fig. 7C) flash rates. Good model fits (dark solid lines) are observed in both SC and V1 at all stimulus rates. Table 2 lists the estimated component amplitudes for both regions and all flash rates (data from 6 animals).

Each phase of the V1 response shows strong dependence on flash rate, whereas the SC response is generally less frequency dependent. In V1, the onset and sustained response components diminish as a function of increasing frequency, whereas the opposite is true for the offset transient. The sustained response attenuates to the extreme when the flash rate reaches 10 Hz leading to potential misinterpretation of the underlying neural responses. By contrast, the SC responses are dominated by sustained signal at all frequencies, with smaller contributions from the transients. A small onset transient that decreases in amplitude with higher flash rates is observed in some SC traces. The offset transient in SC seems to be moderately increasing with flash rate, but unlike in V1, it does not reach in amplitude the signal level reached during sustained stimulation. The poststimulus signal level elevation is only observed in SC at the higher flash rates.

## Discussion

Neurovascular coupling is commonly modeled as a linear time-invariant system that is fully characterized by its IRF (Boynton et al. 1996; Engel et al. 1997). GLM of the BOLD response relies on this assumption and can only accommodate certain special types of nonlinearities with uncertain physiological relevance (Friston et al. 1998). Because interpreting fMRI data in an acceptable manner requires that we know how stable neurovascular coupling is for varying stimuli and across brain regions, using the rat visual system as a model, we evaluated temporal and spatial variances of IRFs and assessed their implications, if any, on GLM of BOLD responses.

### Summary of Present Findings

We built a visual stimulus system to allow change of light color and/or flash rate of full-field binocular stimuli. Our recording setups for fMRI (with multiple slices) and electrophysiology (with multiple electrodes) enabled reproducible measurements of evoked responses from multiple loci (i.e., V1 and SC). Both in SC and V1, reliable and reproducible single-trial responses were obtained with green and blue colors, whereas red color induced only weaker responses. SC time courses showed a rapid rise to peak and moderate decrease during sustained stimulation, whereas V1 exhibited onset/offset transients in addition to the sustained signal during the stimulus. These salient features were depicted in both neural and BOLD responses, albeit at different time scales. Overall, we found that BOLD and neural responses of V1 and SC, which are color sensitive, differ greatly in their temporal evolutions and dependencies on flash rates (Figs 1–3). While IRFV1 and IRFSC were found to be similar at the 1-Hz flash rate, we determined that in the SC temporal invariance generally holds for the flash rate range examined (1–10 Hz), whereas this assumption breaks down significantly in V1 at 10 Hz mainly due to strong attenuation of the sustained BOLD component (i.e., where MUA > BOLD). However, there is a small uncoupling in SC at 10-Hz flash rate for the offset response (i.e., where BOLD > MUA). Both IRFV1 and IRFSC are space invariant within experimental accuracy (Figs 4–6). We illustrate how the use of a suitably designed GLM will accurately represent the hemodynamics at all stimulus rates in both regions, but will lead to the wrong conclusion regarding the underlying neural activity (Fig. 7).

### Comparison with Previous Imaging Studies of the Rat Visual System

Neuroimaging of the rat visual system has recently been performed by optical imaging and fMRI techniques. Gias et al. (2005) demonstrated retinotopic organization in upper layers of rat V1 with optical imaging. To reach deeper cortical layers in the rat, van Camp et al. (2006) used 7.0T fMRI to image functional responses in V1/V2, SC, and parts of the accessory optic system of pigmented Long-Evans rats under medetomidine anesthesia. Although the relatively weak BOLD responses necessitated long stimulus blocks (2 min) and low temporal resolution (25 s), the most salient feature of their results was the opposite frequency dependence of the sustained V1 versus SC responses, a finding which the present study confirms. Furthermore, the authors point to the observation that at high flash rates the SC responses exceed the stimulation period, which is in line with our own observations.

At 9.4T, Pawela et al. (2008) were able to considerably improve the temporal resolution (2 s) and, presumably using a much larger RF coil, covered extensively the posterior brain regions and thereby detected BOLD responses from DLG, V1/V2, and SC in albino Sprague-Dawley rats under medetomidine anesthesia. Signal averaging was necessary across trials to reveal transient responses (i.e., at onset/offset of stimuli) in V1/V2 compared with SC/DLG, a finding which we are also able to corroborate. We note that the modeling approach taken by Pawela et al., which attempts to explain measured BOLD signals as a function of assumed underlying neural dynamics (including refractory nonlinearities), exhibits poorest performance in visual cortex at the higher (5 and 10 Hz) flash rates. We extend this finding by showing that measured neural activity becomes increasingly decoupled from the BOLD response as a function of increasing flash rate, the greatest extent of which was observed at 10-Hz flash rate.

In a recent paper, Lau et al. report on BOLD modulations by short (1-s duration) blocks of 10-Hz flash rate stimulation in SC and thalamus of albino Sprague-Dawley rats under isoflurane anesthesia (Lau et al. 2011). The IRFSC estimated in the present work for 30 s blocks of 1-Hz flash rate stimulation predicts well the SC responses at 5- and 10-Hz flash rates as well. To compare our results with those of Lau et al., we generated BOLD responses by applying IRFSC to the SC MUA data during the first second of 10 Hz stimulation (data not shown). The modeled BOLD response reached half-maximum (“t50”) at 2.8 (±0.2) s and peaked at 4.1 (±0.2) s, which are in excellent agreement with the results obtained by Lau et al. (t50 = 2.6 ± 0.1 s, peak = ∼4 s). The apparently weaker V1 BOLD responses in the Lau et al. study is also in line with our own observations of reduced initial V1 BOLD transients to 10 Hz flash rate stimulation.

The similarities between the works cited above and the present results are extensive, despite the use of varying light delivery strategies (stroboscopic and LED based), varying experimental settings (dark adapted and nonadapted), varying anesthetics (medetomidine, isoflurane, and urethane) and different rat strains (Sprague-Dawley and Long-Evans). Pigmented rats, such as the Long-Evans strain, exhibit some neuroanatomical and functional differences compared with albino strains such as Sprague-Dawley. Albino rats have impaired vision as compared with normally pigmented rats (melanin prevents reflections within the eyeball), and a translucent iris, which leads to retinal degeneration. Albino rats take much longer than pigmented rats to adapt to low-light conditions, have much poorer visual acuity of (20/1200) than pigmented rats (20/600), and are thus, more visually impaired compared with pigmented rats (Lund et al. 1974; Dräger and Olsen 1980; Lund et al. 1980; Prusky et al. 2002). We chose urethane to provide a stable state (Kannurpatti and Biswal 2006) during which sensory-evoked signals are unaffected (Maggi and Meli 1986; Sceniak and Maciver 2006), thus pointing to the compound's suitability in fMRI studies (Yang et al. 1998; Huttunen et al. 2008; Huttunen et al. 2011). Medetomidine, an α2-adrenoreceptor agonist, is thought to produce a lighter state of sedation (Sanganahalli et al. 2009a, 2010). Thus slight dissimilarities in stimulus-evoked magnitudes could partly be due to differences in the baseline excitability levels attained with the 2 anesthetics used, for example, see Maandag et al. (2007).

### Color Vision and Acuity in Rats

The flexible stimulus presentation system allowed us to examine single-trial color dependence of BOLD and neural responses in SC and V1. The rat retina contains 3 classes of photoreceptors with differing sensitivities to light color: rod sensitivity peaks at ∼500 nm (green/cyan), short-wavelength cones peak at ∼360 nm (ultraviolet), and medium-wavelength cones peak at ∼510 nm (green) (Jacobs et al. 2001). Colors used in the present study were blue (455 nm), green (530 nm), and red (627 nm). The BOLD and neural responses to blue and green stimuli were robust, whereas red elicited weaker responses, presumably because the sensitivity peak of the rods and cones are far from the applied red wavelength. Although it remains unclear what functional advantage color vision might impart on a nocturnal species such as the rat, it was shown using behavioral tasks that after training rats were able to discriminate stimuli solely on the basis of their wavelengths, that is, without brightness cues (Jacobs et al. 2001). Recent work by La Garza et al. demonstrate the use of fMRI in the study of the retina, but further work is needed to characterize the hemodynamic and electrophysiological consequences of variations in stimulus properties on downstream processing structures (La Garza et al. 2011).

### Neurovascular Uncouplings in V1 and SC and Implications for GLM

The metabolic and hemodynamic properties of V1 and SC differ (Iadecola et al. 1983; Klein et al. 1986; Maeda et al. 1998). Although the capillary density is lower in SC than in V1, blood flow and energy metabolism in SC are both higher. Studies of cortical single-cell receptive field properties in both rat (Girman et al. 1999; Girman and Lund 2007) and mouse (Niell and Stryker 2008) have shown that—although less organized than those of highly visual mammals (Ohki et al. 2005)—rodent V1 and SC cells are well tuned to orientation and hence offer a well-characterized test bed for neuroscientific inquiry (Binns 1999; Sefton et al. 2004). Afferent signals to V1 are relayed by the thalamus (DLG), whereas SC receives direct retinal innervations. The cortex and thalamus form a reciprocally connected network in which each structure affects the other (Castro-Alamancos 2004; Sherman 2005). Rapid adaptation in V1 after a short initial transient may reflect relatively stronger inhibitory feedback mechanisms (Iadecola 2004), whether they are local or distant, in the cortical networks than in SC.

The majority of the dynamic phases observed in the fMRI data from SC and V1 were confirmed by similar temporal features in the electrophysiology data. Generally, these complimentary observations of BOLD and neural responses indicate linear and time-invariant neurovascular coupling in SC and V1 for most stimuli examined. Furthermore, our results indicate that the couplings in the 2 regions are similar. The sustained V1 and offset SC BOLD responses, however, were uncoupled from the neural activity at the 10-Hz flash rate (i.e., MUA > BOLD in V1 vs. BOLD > MUA in SC). This mismatch between BOLD and neural responses reveals the possibility of neurovascular uncoupling at high flash rates (more so in V1 than in SC), the underlying mechanisms for which could be metabolic, hemodynamic, and/or neurophysiologic in nature. Because we do not fully understand the mechanism for the absence of sustained BOLD signal in V1 despite strong underlying elevated neural activity (and perhaps also the presence of the offset BOLD signal in SC despite small neuronal activity), and both observations which have not yet been observed in humans, future fMRI and electrophysiology studies in conjunction with localized pharmacological intervention in V1, DLG, and/or SC should address these issues.

The use of piecewise linear models in the framework of GLM have been advocated for the study of auditory steady-state responses in humans to account for observed nonlinear responses (Harms and Melcher 2003). To date, there are less clear accounts of similar approaches for the visual modality. Based on the initial observations of the data, we aimed to ascertain whether a GLM consisting of 4 components was able to accurately explain observed dynamics using a single HRF. Our findings indicate that although the GLM provides the desired fit to the data, the results do not represent the underlying neural activity. The spatial distributions of the onset/offset/sustained response components of the cortex strongly overlap, thereby suggesting that most probably the same V1 (and V2) neural networks gives rise to all of these BOLD response components.

A recent human auditory fMRI study showed that high-frequency trains of click sounds produced weaker BOLD signals than predicted by magnetic measurements of the sustained electrical events in primary auditory cortex (Gutschalk et al. 2010). It is important to consider that neither intracortical measures of electrical network activity (in present study) nor extracranial recordings (in the Gutschalk study) can unambiguously distinguish between membrane fluctuations of excitatory versus inhibitory cells. Rather, the power recorded by a microelectrode (in present study) or magnetic sensor (in the Gutschalk study) is a weighted sum of many subpopulations of neurons and possibly even glia. Thus, it remains to be shown whether the weaker BOLD response in V1 at high flash rates is due to a large metabolic demand or a small blood flow response. Preliminary results using surface laser-Doppler flowmetry probes over V1 indicate that blood flow time courses closely mimic the BOLD response at all flash rates examined (data now shown), thereby suggesting that it might be blood flow that is uncoupled from neural activity during the sustained response at high flash rates. While these preliminary results are not as comprehensive as our BOLD and electrophysiological data, we believe that with future improvements in fMRI sensitivity, we can conduct reliable “calibrated fMRI” studies (Herman et al. 2009; Sanganahalli et al. 2009b) of the visual system even though the typical BOLD responses in this sensory system is several times smaller than typical BOLD responses in other sensory systems (Sanganahalli et al. 2008; Sanganahalli, Bailey, et al. 2009).

### Heterogeneity of BOLD Transfer Functions

Several studies have used system identification techniques on human fMRI data to characterize temporal and spatial variability of HRFs (Boynton et al. 1996; Engel et al. 1997; Vazquez and Noll 1998; Birn et al. 2001). Though insightful on stimulus-BOLD coupling, they should be augmented by direct neurovascular coupling analysis using IRFs, which relates local changes in neural activity and dilation (or constriction) of blood vessels to increase (or decrease) blood flow (or volume) to a region. Many animal studies have measured both the neural and vascular responses and show linear neurovascular coupling within some range of stimuli (Logothetis et al. 2001; Matsuura and Kanno 2001; Nemoto et al. 2004; Sheth et al. 2004; Ureshi et al. 2004; Martin et al. 2006; Herman et al. 2009; Sanganahalli et al. 2009b). For example, time invariance of the BOLD transfer function in the somatosensory cortex has been demonstrated (Herman et al. 2009; Sanganahalli et al. 2009b; Hirano et al. 2011) by a procedure based on stimuli of the same amplitude but varying frequencies, as done in the present. While a detailed discussion of each transfer function is not within the present scope, the general observation is that the BOLD transfer function of the somatosensory cortex is generally similar to those obtained here for V1 and SC (Herman et al. 2009; Sanganahalli et al. 2009b).

Much less examined in fMRI is the spatial relationship between neural activity and BOLD signal, which has been measured mainly in the V1 and whisker barrel field of small animals (Kim et al. 2004; Berwick et al. 2008). While these animal studies generally show reasonable neurovascular coupling across small distances in homologous cortical regions, other human studies (which usually do not have measured neural activity) show subtle variations of neurovascular coupling across larger distances in nonhomologous cortical regions (Engel et al. 1997; Birn et al. 2001; Janz et al. 2001; Liu and He 2008; Zhang et al. 2008; Liu et al. 2010; Ostwald et al. 2010; Carlson et al. 2011). However, it remains unclear if these variations have practical consequences for BOLD signal modeling, as differences in experimental noise across the regions have not been fully assessed. In the present study, we found that BOLD transfer functions for SC and V1 are interchangeable. The implication for GLM of BOLD responses is that the use of a single canonical HRF for whole-brain evaluations is likely to be accurate within experimental error, once temporal dynamics of the responses are sufficiently accounted for by the model structure.

## Conclusions

Interpretation of fMRI data mandates that we know how stable neurovascular coupling is for varying stimuli and across brain regions (Boynton et al. 1996; Engel et al. 1997). This is achieved by experimentally measuring IRFs that are calculated from measured neural and BOLD responses. The underlying assumption of linearity and time invariance should not be taken for granted, especially in human studies without neural recordings, though our results suggest that the use of a space-invariant HRF may be sufficient in GLM. Future efforts using dynamic calibrated fMRI (Herman et al. 2009; Sanganahalli et al. 2009b) should be directed to specifically assess neurovascular as well as neurometabolic couplings in a layer-specific manner so that the neurophysiological basis of fMRI may be improved to relate to classical metabolic-flow coupling studies with autoradiographic techniques (Ueki et al. 1988, 1992).

## Supplementary Material

Supplementary material can be found at: http://www.cercor.oxfordjournals.org/

## Funding

The work was supported by grants from National Institutes of Health (R01 MH067528 and P30 NS52519 to F.H.; R01 NS-049307 and R01 NS-066974 to H.B.). C.J.B. is grateful to Aarhus University for his graduate program scholarship and to Dagmar Marshall's Fond, Oticon Fonden, Brorsons Rejselegat, and Torben and Alice Frimodts Fond for financial support.

The authors thank colleagues (Peter Brown, Scott McIntyre, Bei Wang, and Alyssa Siefert) at MRRC (mrrc.yale.edu) and Core Center for QNMR (qnmr.yale.edu) for assistance, materials, and help with building the experimental setup. Conflict of Interest: None declared.

## References

Berwick
J
Johnston
D
Jones
M
Martindale
J
Martin
C
Kennerley
AJ
Redgrave
P
Mayhew
JE
Fine detail of neurovascular coupling revealed by spatiotemporal analysis of the hemodynamic response to single whisker stimulation in rat barrel cortex
J Neurophysiol
,
2008
, vol.
99
(pg.
787
-
798
)
Binns
KE
The synaptic pharmacology underlying sensory processing in the superior colliculus
Prog Neurobiol
,
1999
, vol.
59
(pg.
129
-
159
)
Birn
RM
ZS
Bandettini
PA
Spatial heterogeneity of the nonlinear dynamics in the FMRI BOLD response
Neuroimage
,
2001
, vol.
14
(pg.
817
-
826
)
Borowsky
IW
Collins
RC
Metabolic anatomy of brain: a comparison of regional capillary density, glucose metabolism, and enzyme activities
J Comp Neurol
,
1989
, vol.
288
(pg.
401
-
413
)
Boynton
GM
Engel
SA
Glover
GH
Heeger
DJ
Linear systems analysis of functional magnetic resonance imaging in human V1
J Neurosci
,
1996
, vol.
16
(pg.
4207
-
4221
)
Carlson
T
Hogendoorn
H
Fonteijn
H
Verstraten
FA
Spatial coding and invariance in object-selective cortex
Cortex
,
2011
, vol.
47
(pg.
14
-
22
)
Castro-Alamancos
MA
Dynamics of sensory thalamocortical synaptic networks during information processing states
Prog Neurobiol
,
2004
, vol.
74
(pg.
213
-
247
)
Dräger
UC
Olsen
JF
Origins of crossed and uncrossed retinal projections in pigmented and albino mice
J Comp Neurol
,
1980
, vol.
191
(pg.
383
-
412
)
Engel
SA
Glover
GH
Wandell
BA
Retinotopic organization in human visual cortex and the spatial precision of functional MRI
Cereb Cortex
,
1997
, vol.
7
(pg.
181
-
192
)
Englot
DJ
Mishra
AM
Mansuripur
PK
Herman
P
Hyder
F
Blumenfeld
H
Remote effects of focal hippocampal seizures on the rat neocortex
J Neurosci
,
2008
, vol.
28
(pg.
9066
-
9081
)
Frackowiak
RSJ
Human brain function
,
2004

Frahm
J
Haase
A
Matthaei
D
Rapid three-dimensional MR imaging using the FLASH technique
J Comput Assist Tomogr
,
1986
, vol.
10
(pg.
363
-
368
)
Friston
KJ
Imaging neuroscience: principles or maps?
Proc Natl Acad Sci U S A
,
1998
, vol.
95
(pg.
796
-
802
)
Friston
KJ
Josephs
O
Rees
G
Turner
R
Nonlinear event-related responses in fMRI
Magn Reson Med
,
1998
, vol.
39
(pg.
41
-
52
)
Gias
C
Hewson-Stoate
N
Jones
M
Johnston
D
Mayhew
JE
Coffey
PJ
Retinotopy within rat primary visual cortex using optical imaging
Neuroimage
,
2005
, vol.
24
(pg.
200
-
206
)
Girman
SV
Lund
RD
Most superficial sublamina of rat superior colliculus: neuronal response properties and correlates with perceptual figure-ground segregation
J Neurophysiol
,
2007
, vol.
98
(pg.
161
-
177
)
Girman
SV
Sauvé
Y
Lund
RD
Receptive field properties of single neurons in rat primary visual cortex
J Neurophysiol
,
1999
, vol.
82
(pg.
301
-
311
)
Goense
JBM
Logothetis
NK
Neurophysiology of the BOLD fMRI signal in awake monkeys
Curr Biol
,
2008
, vol.
18
(pg.
631
-
640
)
Gutschalk
A
Hamalainen
MS
Melcher
JR
BOLD responses in human auditory cortex are more closely related to transient MEG responses than sustained ones
J Neurophysiol
,
2010
, vol.
103
(pg.
2015
-
2026
)
Harms
MP
Melcher
JR
Detection and quantification of a wide range of fMRI temporal responses using a physiologically-motivated basis set
Hum Brain Mapp
,
2003
, vol.
20
(pg.
168
-
183
)
Hennig
J
Nauerth
A
Friedburg
H
RARE imaging: a fast imaging method for clinical MR
Magn Reson Med
,
1986
, vol.
3
(pg.
823
-
833
)
Herman
P
Sanganahalli
BG
Blumenfeld
H
Hyder
F
Cerebral oxygen demand for short-lived and steady-state events
J Neurochem
,
2009
, vol.
109

Suppl 1
(pg.
73
-
79
)
Hespanha
JP
Linear systems theory
,
2009
Princeton (NJ)
Princeton University Press
Hirano
Y
Stefanovic
B
Silva
AC
Spatiotemporal evolution of the functional magnetic resonance imaging response to ultrashort stimuli
J Neurosci
,
2011
, vol.
31
(pg.
1440
-
1447
)
Huettel
SA
Song
AW
McCarthy
G
Functional magnetic resonance imaging
,
2004
Sunderland (MA)
Sinauer Associates, Inc
Huttunen
JK
Gröhn
O
Penttonen
M
Coupling between simultaneously recorded BOLD response and neuronal activity in the rat somatosensory cortex
Neuroimage
,
2008
, vol.
39
(pg.
775
-
785
)
Huttunen
JK
Niskanen
J-P
Lehto
LJ
Airaksinen
AM
Niskanen
EI
Penttonen
M
Gröhn
O
Evoked local field potentials can explain temporal variation in blood oxygenation level-dependent responses in rat somatosensory cortex
NMR Biomed
,
2011
, vol.
24
(pg.
209
-
215
)
Hyder
F
Kida
I
Behar
KL
Kennan
RP
Maciejewski
PK
Rothman
DL
Quantitative functional imaging of the brain: towards mapping neuronal activity by BOLD fMRI
NMR Biomed
,
2001
, vol.
14
(pg.
413
-
431
)
Hyder
F
Rothman
DL
Blamire
AM
Image reconstruction of sequentially sampled echo-planar data
Magn Reson Imaging
,
1995
, vol.
13
(pg.
97
-
103
)
Hyder
F
Sanganahalli
BG
Herman
P
Coman
D
Maandag
NJG
Behar
KL
Blumenfeld
H
Rothman
DL
Neurovascular and neurometabolic couplings in dynamic calibrated fMRI: transient oxidative neuroenergetics for block-design and event-related paradigms
Front Neuroenergetics
,
2010
, vol.
2

C
Neurovascular regulation in the normal brain and in Alzheimer's disease
Nat Rev Neurosci
,
2004
, vol.
5
(pg.
347
-
360
)
C
Nakai
M
Mraovitch
S
Ruggiero
DA
Tucker
LW
Reis
DJ
Global increase in cerebral metabolism and blood flow produced by focal electrical stimulation of dorsal medullary reticular formation in rat
Brain Res
,
1983
, vol.
272
(pg.
101
-
114
)
Jacobs
GH
Fenwick
JA
Williams
GA
Cone-based vision of rats for ultraviolet and visible lights
J Exp Biol
,
2001
, vol.
204
(pg.
2439
-
2446
)
Janz
C
Heinrich
SP
Kornmayer
J
Bach
M
Hennig
J
Coupling of neural activity and BOLD fMRI response: new insights by combination of fMRI and VEP experiments in transition from single events to continuous stimulation
Magn Reson Med
,
2001
, vol.
46
(pg.
482
-
486
)
Kannurpatti
SS
Biswal
BB
Spatial extent of CBF response during whisker stimulation using trial averaged laser Doppler imaging
Brain Res
,
2006
, vol.
1089
(pg.
135
-
142
)
Kim
D-S
Ronen
I
Olman
C
Kim
S-G
Ugurbil
K
Toth
LJ
Spatial relationship between neuronal activity and BOLD functional MRI
Neuroimage
,
2004
, vol.
21
(pg.
876
-
885
)
Klein
B
Kuschinsky
W
Schrock
H
Vetterlein
F
Interdependency of local capillary density, blood flow, and metabolism in rat brains
Am J Physiol
,
1986
, vol.
251
(pg.
H1333
-
H1340
)
La Garza
BHD
Muir
ER
Li
G
Shih
Y-YI
Duong
TQ
Blood oxygenation level-dependent (BOLD) functional MRI of visual stimulation in the rat retina at 11.7 T
NMR Biomed
,
2011
, vol.
24
(pg.
183
-
193
)
Lau
C
Zhou
IY
Cheung
MM
Chan
KC
Wu
EX
BOLD temporal dynamics of rat superior colliculus and lateral geniculate nucleus following short duration visual stimulation
PLoS One
,
2011
, vol.
6
pg.
e18914

Liu
Z
He
B
fMRI-EEG integrated cortical source imaging by use of time-variant spatial constraints
Neuroimage
,
2008
, vol.
39
(pg.
1198
-
1214
)
Liu
Z
Rios
C
Zhang
N
Yang
L
Chen
W
He
B
Linear and nonlinear relationships between visual stimuli, EEG and BOLD fMRI signals
Neuroimage
,
2010
, vol.
50
(pg.
1054
-
1066
)
Liu
Z
Zhang
N
Chen
W
He
B
Mapping the bilateral visual integration by EEG and fMRI
Neuroimage
,
2009
, vol.
46
(pg.
989
-
997
)
Logothetis
NK
Pauls
J
Augath
M
Trinath
T
Oeltermann
A
Neurophysiological investigation of the basis of the fMRI signal
Nature
,
2001
, vol.
412
(pg.
150
-
157
)
Lund
RD
Land
PW
Boles
J
Normal and abnormal uncrossed retinotectal pathways in rats: an HRP study in adults
J Comp Neurol
,
1980
, vol.
189
(pg.
711
-
720
)
Lund
RD
Lund
JS
Wise
RP
The organization of the retinal projection to the dorsal lateral geniculate nucleus in pigmented and albino rats
J Comp Neurol
,
1974
, vol.
158
(pg.
383
-
403
)
Maandag
NJG
Coman
D
Sanganahalli
BG
Herman
P
Smith
AJ
Blumenfeld
H
Shulman
RG
Hyder
F
Energetics of neuronal signaling and fMRI activity
Proc Natl Acad Sci U S A
,
2007
, vol.
104
(pg.
20546
-
20551
)
M
A simplified formulation of the gamma variate function
Phys Med Biol
,
1992
, vol.
37
(pg.
1597
-
1600
)
Maeda
M
Duelli
R
Schrock
H
Kuschinsky
W
Autoradiographic determination of local cerebral blood flow and local cerebral glucose utilization during chemical stimulation of the nucleus tractus solitarii of anesthetized rats
J Auton Nerv Syst
,
1998
, vol.
69
(pg.
132
-
140
)
Maggi
CA
Meli
A
Suitability of urethane anesthesia for physiopharmacological investigations in various systems. Part 1: general considerations
Experientia
,
1986
, vol.
42
(pg.
109
-
114
)
Martin
C
Martindale
J
Berwick
J
Mayhew
J
Investigating neural-hemodynamic coupling and the hemodynamic response function in the awake rat
Neuroimage
,
2006
, vol.
32
(pg.
33
-
48
)
Matsuura
T
Kanno
I
Quantitative and temporal relationship between local cerebral blood flow and neuronal activation induced by somatosensory stimulation in rats
Neurosci Res
,
2001
, vol.
40
(pg.
281
-
290
)
Nemoto
M
Sheth
S
Guiou
M
Pouratian
N
Chen
JW
Toga
AW
Functional signal- and paradigm-dependent linear relationships between synaptic activity and hemodynamic responses in rat somatosensory cortex
J Neurosci
,
2004
, vol.
24
(pg.
3850
-
3861
)
Niell
CM
Stryker
MP
Highly selective receptive fields in mouse visual cortex
J Neurosci
,
2008
, vol.
28
(pg.
7520
-
7536
)
Ogawa
S
Lee
TM
Kay
AR
Tank
DW
Brain magnetic resonance imaging with contrast dependent on blood oxygenation
Proc Natl Acad Sci U S A
,
1990
, vol.
87
(pg.
9868
-
9872
)
Ohki
K
Chung
S
Ch'ng
YH
Kara
P
Reid
RC
Functional imaging with cellular resolution reveals precise micro-architecture in visual cortex
Nature
,
2005
, vol.
433
(pg.
597
-
603
)
Ostwald
D
Porcaro
C
Bagshaw
AP
An information theoretic approach to EEG-fMRI integration of visually evoked responses
Neuroimage
,
2010
, vol.
49
(pg.
498
-
516
)
Pawela
CP
Hudetz
AG
Ward
BD
Schulte
ML
Li
R
Kao
DS
Mauck
MC
Cho
YR
Neitz
J
Hyde
JS
Modeling of region-specific fMRI BOLD neurovascular response functions in rat brain reveals residual differences that correlate with the differences in regional evoked potentials
Neuroimage
,
2008
, vol.
41
(pg.
525
-
534
)
Paxinos
G
Watson
C
The rat brain in stereotaxic coordinates
,
1998
San Diego (CA)
Prusky
GT
Harker
KT
Douglas
RM
Whishaw
IQ
Variation in visual acuity within pigmented, and between pigmented and albino rat strains
Behav Brain Res
,
2002
, vol.
136
(pg.
339
-
348
)
Sanganahalli
BG
Herman
P
Blumenfeld
H
Hyder
F
fMRI and electrophysiological studies with α-chloralose and domitor anesthesia
J Cereb Blood Flow Metab
,
2009
, vol.
29
pg.
S616

Sanganahalli
BG
Herman
P
Blumenfeld
H
Hyder
F
BOLD impulse response functions and baseline-dependent response adaptation
Proc Int Soc Magn Reson Med
,
2010
, vol.
18
pg.
184

Sanganahalli
BG
Bailey
CJ
Herman
P
Hyder
F
Tactile and non-tactile sensory paradigms for fMRI and neurophysiologic studies in rodents
Methods Mol Biol
,
2009
, vol.
489
(pg.
213
-
242
)
Sanganahalli
BG
Herman
P
Blumenfeld
H
Hyder
F
J Neurosci
,
2009
, vol.
29
(pg.
1707
-
1718
)
Sanganahalli
BG
Herman
P
Hyder
F
Frequency-dependent tactile responses in rat brain measured by functional MRI
NMR Biomed
,
2008
, vol.
21
(pg.
410
-
416
)
Sceniak
MP
Maciver
MB
Cellular actions of urethane on rat visual cortical neurons in vitro
J Neurophysiol
,
2006
, vol.
95
(pg.
3865
-
3874
)
Sefton
AJ
Dreher
B
Harvey
A
Paxinos
G
Visual system
The rat nervous system
,
2004
Amsterdam
(pg.
1083
-
1165
)
Sherman
SM
Thalamic relays and cortical functioning
Prog Brain Res
,
2005
, vol.
149
(pg.
107
-
126
)
Sheth
SA
Nemoto
M
Guiou
M
Walker
M
Pouratian
N
Toga
AW
Linear and nonlinear relationships between neuronal activity, oxygen metabolism, and hemodynamic responses
Neuron
,
2004
, vol.
42
(pg.
347
-
355
)
Smith
AJ
Blumenfeld
H
Behar
KL
Rothman
DL
Shulman
RG
Hyder
F
Cerebral energetics and spiking frequency: the neurophysiological basis of fMRI
Proc Natl Acad Sci U S A
,
2002
, vol.
99
(pg.
10765
-
10770
)
Stark
E
Abeles
M
Predicting movement from multiunit activity
J Neurosci
,
2007
, vol.
27
(pg.
8387
-
8394
)
Ueki
M
Linn
F
Hossmann
KA
Functional activation of cerebral blood flow and metabolism before and after global ischemia of rat brain
J Cereb Blood Flow Metab
,
1988
, vol.
8
(pg.
486
-
494
)
Ueki
M
Mies
G
Hossmann
KA
Effect of alpha-chloralose, halothane, pentobarbital and nitrous oxide anesthesia on metabolic coupling in somatosensory cortex of rat
Acta Anaesthesiol Scand
,
1992
, vol.
36
(pg.
318
-
322
)
Ureshi
M
Matsuura
T
Kanno
I
Stimulus frequency dependence of the linear relationship between local cerebral blood flow and field potential evoked by activation of rat somatosensory cortex
Neurosci Res
,
2004
, vol.
48
(pg.
147
-
153
)
van Camp
N
Verhoye
M
de Zeeuw
CI
van der Linden
A
Light stimulus frequency dependence of activity in the rat visual system as studied with high-resolution BOLD fMRI
J Neurophysiol
,
2006
, vol.
95
(pg.
3164
-
3170
)
van der Linden
A
van Camp
N
Ramos-Cabrer
P
Hoehn
M
Current status of functional MRI on small animals: application to physiology, pathophysiology, and cognition
NMR Biomed
,
2007
, vol.
20
(pg.
522
-
545
)
Vazquez
AL
Noll
DC
Nonlinear aspects of the BOLD response in functional MRI
Neuroimage
,
1998
, vol.
7
(pg.
108
-
118
)
Yang
X
Renken
R
Hyder
F
Siddeek
M
Greer
CA
Shepherd
GM
Shulman
RG
Dynamic mapping at the laminar level of odor-elicited responses in rat olfactory bulb by functional MRI
Proc Natl Acad Sci U S A
,
1998
, vol.
95
(pg.
7715
-
7720
)
Zhang
N
Liu
Z
He
B
Chen
W
Noninvasive study of neurovascular coupling during graded neuronal suppression
J Cereb Blood Flow Metab
,
2008
, vol.
28
(pg.
280
-
290
)