Long-Latency Reductions in Gamma Power Predict Hemodynamic Changes That Underlie the Negative BOLD Signal

Studies that use prolonged periods of sensory stimulation report associations between regional reductions in neural activity and negative blood oxygenation level-dependent (BOLD) signaling. However, the neural generators of the negative BOLD response remain to be characterized. Here, we use single-impulse electrical stimulation of the whisker pad in the anesthetized rat to identify components of the neural response that are related to “negative” hemodynamic changes in the brain. Laminar multiunit activity and local field potential recordings of neural activity were performed concurrently with two-dimensional optical imaging spectroscopy measuring hemodynamic changes. Repeated measurements over multiple stimulation trials revealed significant variations in neural responses across session and animal datasets. Within this variation, we found robust long-latency decreases (300 and 2000 ms after stimulus presentation) in gamma-band power (30–80 Hz) in the middle-superficial cortical layers in regions surrounding the activated whisker barrel cortex. This reduction in gamma frequency activity was associated with corresponding decreases in the hemodynamic responses that drive the negative BOLD signal. These findings suggest a close relationship between BOLD responses and neural events that operate over time scales that outlast the initiating sensory stimulus, and provide important insights into the neurophysiological basis of negative neuroimaging signals.


Introduction
Blood oxygenation level-dependent (BOLD) functional magnetic resonance imaging (fMRI) has revolutionized our ability to investigate human brain function noninvasively. The BOLD signal, however, provides only an indirect measure of neuronal activity. This measured secondary hemodynamic marker is linked to associated neural responses via neurovascular coupling. Positive BOLD signals are robustly coupled to local increases in synaptic activity in the form of local field potentials (LFPs) and multiunit activity (MUA;Logothetis, 2008). In contrast, negative BOLD signals are poorly understood (Lauritzen et al., 2012;Moraschi et al., 2012). Although sparse, recent evidence has shown negative BOLD signals to be associated with local decreases in LFP and MUA (Shmuel et al., 2006;Boorman et al., 2010), sug-gesting a neurogenic, rather than vasogenic, basis of generation (Shmuel et al., 2006;Devor et al., 2007;Boorman et al., 2010;Liu et al., 2011;Yuan et al., 2011;Mullinger et al., 2014). However, the component of the neural response that is responsible for negative perfusion response-related signals remains uncharacterized and is essential to studying the balance between neural excitation and inhibition in and across different brain regions using neuroimaging tools. We therefore sought to identify the components in the neural response that underpin negative hemodynamics using a multimodal neuroimaging approach (Logothetis, 2008;Liu et al., 2011;Moraschi et al., 2012). Using our recently established model of negative BOLD in the anesthetized rodent , along with concurrent two-dimensional optical imaging spectroscopy (2D-OIS) and multichannel/multiprobe electrophysiology, we examined the neurophysiological basis of positive and negative hemodynamic responses evoked in primary somatosensory cortex by electrical stimulation of the whisker pad. Two-dimensional OIS provides high-resolution measures of changes in total hemoglobin and saturation, which underpin the BOLD signal (Kennerley et al., 2012b). Previous research by our laboratory and those of others have typically used long-duration stimulation paradigms, such as flashing checker boards (Shmuel et al., 2006) or trains of stimulation  to investigate the negative BOLD signal. However, while such paradigms elicit robust and easily measurable signal changes, repeated stimuli generate hemodynamic responses that cannot be linearly decomposed using impulse response functions obtained from single-impulse stimuli , rendering interpretation of the data and underlying neurovascular relationships nontrivial. Therefore, in the current study, we presented singleimpulse stimulation alongside established long-duration "train" stimuli to examine the value of each form of stimulation when identifying neural markers of negative hemodynamic signals. Our principal finding was that the negative hemodynamic response that surrounded the regionally restricted positive hemodynamic response was associated with a long-latency (300 -2000 ms after stimulation) reduction of gamma power in the evoked neural response to single-impulse stimulation. This gamma power reduction was confined to the middle-superficial cortical layers. Interestingly, this was not observed during long-duration stimulation, suggesting a masking effect of repeated stimuli. Our results highlight the value of single-impulse stimulation to investigate neurovascular coupling, provide significant insights into the etiology of negative perfusion response-related imaging signals, and have important implications for the interpretation of BOLD fMRI signals.

Materials and Methods
Experimental design. Multichannel recording electrodes were placed simultaneously into an activated region of whisker barrel cortex and into a surrounding area that exhibited a negative hemodynamic response. Current-regulated electrical stimulation was applied to the whole contralateral whisker pad. Natural variation in evoked neural activity and accompanying hemodynamic responses were used to investigate differences in the neurovascular coupling on a trial-by-trial basis. Individual stimulation trials from each subject were sorted and grouped according to either the magnitude of the evoked raw LFP or the frequency power of the neural activity. The trial-matched concurrently recorded hemodynamic data were then averaged to generate time series of the change in the localized concentration and saturation of hemoglobin. We found that the sensory-evoked surround negative hemodynamic response was associated with an extremely long latency (300 -2000 ms after stimulation) local reduction of gamma power. Implications this has for the interpretation of neuroimaging signals will be discussed.
Animal preparation and surgery. Surgical and experimental procedures were performed in accordance with the Animal (Scientific Procedures) Act 1986, with approval from the United Kingdom Home Office. Female hooded Lister rats (n ϭ 10) were kept at 22°C in a 12 h dark/light cycle with food and water supplied ad libitum. Rats weighing between 230 and 330 g were anesthetized with urethane at 1.25 g/kg, intraperitoneally; additional supplementary doses of 0.1 ml were administered if required. Atropine was administered (0.4 mg/kg, s.c.) to lessen mucous secretions. A homoeothermic heating blanket system with rectal temperature monitoring (Harvard Instruments) maintained core temperature at 37°C throughout surgery and subsequent experimental procedures. Animals were tracheotomized to allow artificial ventilation and continuous monitoring of end-tidal CO 2 (CapStar-100, CWE Systems). Cannulation of femoral arteries was used to monitor mean arterial blood pressure (MABP) and to allow blood samples to be taken for blood gas measurement. The blood gas measurements were used to adjust the ventilator parameters to keep the blood oxygen saturation of the subject within physiological limits. A femoral vein was also cannulated to allow an infusion of phenylephrine (0.13-0.26 mg/h) so that MABP could be maintained between 100 and 110 mmHg (Golanov et al., 1994;Nakai and Maeda, 1999). Subjects were placed in a stereotaxic frame (Kopf Instruments) and the skull surface was exposed. A dental drill, cooled with saline, was used to thin an area of bone overlying the right somatosensory and motor cortices to translucency. Dental cement was used to secure a plastic well (20 mm diameter) with input/output ports over the thinned window. The well was constantly infused with saline (ϳ10 ml/h) to maintain window transparency and to reduce optical specularities reflected from the skull surface.
Electrical stimulation of the whisker pad. Two insulated stainless steel electrodes (exposed tip, 2 mm) were inserted subcutaneously into the left whisker pad between vibrissae rows A/B and C/D. Muscle responses evoked by electrical stimulation (0.8 mA; 300 s pulse width, 2 s duration at 5 Hz) confirmed that most of the whisker pad was stimulated. The electrical stimuli used in this study produced no reliable changes in MABP, CO 2 partial pressure, or heart rate, thereby precluding potentially confounding changes in systemic physiology affecting evoked neural and hemodynamic responses in somatosensory cortex.
Two-dimensional OIS. The region of thinned skull overlying the right somatosensory cortex was illuminated in sequence with four wavelengths of light using a Lambda DG-4 high-speed galvanometer-based filter changer with integrated light source (Sutter Instrument). Images of the cortical surface were recorded by a CCD camera (1M30P, Teledyne DALSA), operating at 4 ϫ 4 binning. Each image pixel represented ϳ75 ϫ 75 m of cortical tissue. The camera had a quantum efficiency of 28% at 500 nm. Image capture was at 32 Hz, and was used to synchronize filter changing for each of the four wavelength exposures, giving an effective frame rate of 8 Hz. The four wavelengths were chosen specifically as two pairs (495 nm Ϯ 31 FWHM and 559 nm Ϯ 16 FWHM; 575 nm Ϯ 14 FWHM, and 587 nm Ϯ 9 FWHM). The wavelengths in each pair have similar total absorption coefficients and therefore sample the same tissue volume. However, they have maximally different absorption coefficients for oxyhemoglobin (HbO 2 ) and deoxyhemoglobin (Hbr) to maximize signal-to-noise ratios. Recorded images were subject to spectral analysis using a path length-scaling algorithm. The algorithm is based on a modified Beer-Lambert law with a path length correction factor . The algorithm requires a baseline value for the concentration of hemoglobin in tissue, which, according to previous measurements , we set to 104 M. Initially, pixel-by-pixel analyses  produced 2D spatiotemporal "maps" of estimated measures of HbO 2 , Hbr, and total hemoglobin (Hbt). Hbt was chosen for subsequent analyses of hemodynamic changes so that, in each trial, the Hbt measurement represents a micromolar change in concentration from its own prestimulus baseline value.
Localization of the whisker barrel and surrounding regions. We used 2D-OIS to provide a functional localization of the whisker barrel cortex by detecting the hemodynamic changes that occurred in response to brief stimulation of the whisker pad [30 stimulation trials; 300 s pulse width; 5 Hz for 2 s; constant current of 0.8 mA; interstimulus interval (ISI) of 25 s]. The intensity of 0.8 mA is in the midrange of stimuli that have been used to elicit cortical neural and hemodynamic responses (Sheth et al., 2004;Franceschini et al., 2008). It is marginally higher than that used to stimulate the whisker pads of awake behaving rodents (0.5 mA; Martin et al., 2006). Images from each stimulus presentation trial were trial averaged before being subjected to spectral analysis. A general linear model (GLM) was applied to 2D-trial-averaged hemodynamic data and the fit of the GLM component representative of the hemodynamic response produced a spatial map of z-score fit. Pixels with the highest z-score and spatially connected with those Ͼ50% of the maximum were included in a "positive" region of interest (Fig. 1F, red). Those with highest anticorrelated z-score fit, again spatially connected and Ͼ50% of the maximum, were included in the "negative" region of interest (Fig. 1F, green). Pixels within the two regions of interest were independently averaged to generate the hemodynamic time series of Hbt, HbO 2 , and Hbr for positive cortical whisker barrel and negative surround regions.
This information was used to guide the placement of two multichannel electrode arrays into the positive whisker barrel somatosensory cortex and into a surround negative region. Care was taken to avoid surface vasculature at the sites of electrode placement. Small holes (diameter, ϳ300 m) were drilled in the thinned skull directly over the points selected for electrode placement. After piercing the dura with a 27 gauge needle, a stereotaxic arm (Kopf Instruments) was used to insert the 16channel silicone linear array electrodes (100 m spacing; area of each site, 177 m 2 ; Neuronexus Technologies) normal to the cortical surface to a depth of 1600 m. The probes were coupled to a preamplifier that was in turn connected to a data acquisition unit via fiber optic cable (Medusa Bioamp, Tucker Davis Technologies). Stimulus-evoked and spontaneously occurring neural recordings were sampled at 24.41 kHz with a 16-bit resolution. Data were collected using OpenEX software (Tucker Davis Technologies), which also handled trigger timings and data storage. Following electrode placement, the whisker pad was again stimulated using the same 2 s train paradigm, now with concurrent recording of neural activity and 2D-OIS. These procedures confirmed both electrode function and appropriate electrode placements within the whisker barrel field and surrounding negative regions.
Experimental paradigms. The whisker pad was stimulated to evoke neural and associated hemodynamic responses. A brief stimulus pulse train of 2 s (5 Hz, 300 s pulse width, 0.8 mA intensity) was presented to localize the regions exhibiting positive and negative hemodynamic responses for subsequent electrode placement (Fig. 1). Three stimulation paradigms were used: one with a long-duration train of stimulating pulses (16 s, 5 Hz, 300 s pulse width, 0.8 mA intensity), one with singlepulse stimulation (300 s pulse width, 0.8 mA intensity), and a recording-only control paradigm where no stimulation was delivered. The order of presenting the four paradigms was randomized for each subject. Neural activity was recorded for 5.2 s before stimulation and for either 6.8 or 21.5 s following the onset of stimulation, depending on the duration of stimulus train while 2D-OIS was recorded throughout each of the stimulus paradigms.
Whisker stimulus presentation paradigms for trial sorting. In animal neuroimaging and human fMRI experiments, it is common practice to record the hemodynamic responses evoked by trains of sensory stimuli. In these experiments, responses to individual stimuli within the train are effectively "summed." This increases the overall magnitude of the hemodynamic response, enhancing the signal-to-noise ratio, and making the response easier to detect. Therefore, in the current study the first stimulus paradigm used stimulus trains (5 Hz) with duration of 16 s and ISI of 70 s (Fig. 2). We showed previously that 30 trials using these conditions produced robust hemodynamic responses measured by BOLD fMRI , 2D-OIS (Kennerley et al., 2012a), and laser Doppler flowmetry (LDF) . In addition, single-impulse stimuli (Hirano et al., 2011) were presented to each subject 100 times with an ISI of 25 s. This condition avoids the complexities of interpreting the responses to multiple neural activations and allows the uninterrupted time course of the neurovascular coupling to be examined, including the complete return to baseline.
Electrophysiological data analysis. Data recorded from the multichannel electrodes were analyzed for both raw LFP activity and activity within specific frequency bands. Recorded LFPs are thought to represent afferent activity and intrinsic neural processing, while MUA represents spiking output (Logothetis, 2008). LFPs represent the superposition of current flowing into and out of activated cells with the summed changes in potential made up from a range of sources. The LFPs are primarily generated from the perisynaptic activity of neurons from a large spatial area, but also include single-unit activity from neurons close to the electrode and MUA spiking activity from local populations of nearby neurons (Logothetis, 2008). Strong correlations have been demonstrated between BOLD signals and LFPs (Logothetis et al., 2001;Goense and Logothetis, 2008), and, while MUA has also been positively correlated with BOLD, BOLD signaling has been observed in the absence of spiking activity (Logothetis, 2008). Other studies have demonstrated further associations between BOLD and specific frequency bands of  Figure 1. Spatial hemodynamic responses evoked by 0.8 mA electrical stimulation of the whisker pad (2 s train of pulses presented at 5 Hz) from four representative animals. A, In vivo grayscale CCD camera image of thinned cranial window, which includes the primary somatosensory and motor cortices. Midline is to the left. B, C, Postmortem histological section (50 m thick) centered on the S1 whisker barrel field and stained for cytochrome oxidase. The whisker barrels and surrounding anatomically mapped regions have been outlined and labeled. The blue overlaid patch represents the approximate location of the negative hemodynamic changes that occurred in response to stimulation of the whole whisker pad. D, In vivo grayscale CCD camera image of thinned cranial window. E, Activation map of change in Hbt, generated by statistical parametric mapping (SPM) GLM with boxcar hemodynamic response function. Scale bar represents z-score. F, Regions of interest selected from activation map in B. Regions having Hbt increases were selected by taking pixels having a z-scores Ͼ50% of the maximum z-score, while regions having Hbt decreases were selected by taking pixels having a z-score Ͻ50% of the minimum z-score. Red represents an increase in Hbt, overlying the activated whisker barrel field. Green represents surrounding regions of decreased Hbt and purple represents the increase in Hbt and overlies the primary motor cortex. G, Outline of multichannel recording electrodes inserted into the whisker barrel region (W) and surround region (S).
Band-limited frequency analysis of LFP recordings offers a basic level of signal disambiguation by selecting only those frequencies within the chosen range. For the present study, classically defined EEG bands were selected for analysis: delta, 1-4 Hz; theta, 4 -8 Hz; alpha, 8 -14 Hz; beta, 14 -30 Hz; and gamma, 30 -80 Hz. As higher frequencies can be observed in extracellular recordings, MUA (300 -3000 Hz) was also examined. Frequency-band analysis was performed using frequency-power analysis (Shmuel et al., 2006;Boorman et al., 2010). Data were placed into 100 ms temporal bins and subjected to Fourier transform as this allowed investigation of the temporal dynamics of the neural activity within the gamma (30 -80 Hz) and MUA (300 -3000 Hz) frequency bands. Once the temporal components of the neural activity had been examined, neural data were taken over time ranges Ͼ100 ms (e.g., 900 and 1700 ms) and subject to Fourier transform, before being separated into the selected range of frequency bands. This allowed frequencies Ͻ10 Hz to be observed, such as in the delta, theta, and alpha bands. For single-trial analysis, data were downsampled ("decimated") to 6.104 kHz for trial sorting based on LFP and MUA and to 1.526 kHz for trial sorting using the classical EEG frequency bands.
Sorting data from individual trials. Trial-by-trial analyses of neural recordings were used to further uncover the underpinnings of the negative neurovascular relationship. This approach relies upon natural variations in the time course and magnitude of the evoked responses on individual trials. These variations are commonplace despite the imposition of identical stimulus conditions. Such variations are frequently deemed detrimental and their unwanted effects on experimental outcome are overcome by averaging across stimulation trials and subjects. In the current study, however, we use these natural response variations to investigate the relationship between negative hemodynamic responses and the underlying changes in neural activity. This approach is commonly used in EEG research, with a number of studies performing trial-by-trial analysis when recording concurrent BOLD and EEG ( Thus, following analysis of electrode recordings, four specific components of the neural response were taken and used to perform trial-by-trial sorting of both the neural activity and their accompanying hemodynamic responses. The components used for sorting were based on either components of the LFP or on the frequency-band power. The LFP component was taken as the initial depolarization that follows each stimulus and was present in the activated region of whisker barrel cortex. The six frequency-band ranges were investigated over specific time periods and spatially across the linearly spaced electrode channels. Spontaneous (nonstimulation) recordings. Trial sorting of individual hemodynamic responses based on concurrently recorded neural activity will be influenced by ongoing spontaneous neural activity, as well as by the neural activity evoked by the sensory stimulation. To determine the effects of spontaneous neural activity on hemodynamic and neural recordings, we included "no-stimulus" trials matched to the 25 s ISI, 100 trial single-impulse stimulation paradigm. This allowed the same analyses to be run on the nonstimulated data, as applied to the stimulation data, to establish the baseline effects of spontaneous neural activation on trial sorting.
Randomized trial sorting. To ensure that the trial-sorting method did not produce any bias by itself generating erroneous differences in the magnitude of sorted hemodynamic responses, additional sorting was performed using randomized vectors of trials. The randomized vectors were applied to both the single-impulse and nonstimulation datasets.
Spatial analysis of sorted hemodynamic data. After the individual trials of hemodynamic data had been sorted according to the various neural criteria, the top 33% and bottom 33% of trials were averaged. This produced a spatial image block over time where regions of interest around each electrode tip (whisker barrel and surround) could be selected and averaged time series generated from these regions. This allowed for temporal analyses of the evoked responses. Spatially averaged hemodynamic data were subject to a concentric ring analysis based on the center of the positively activated whisker region. Circular rings spaced ϳ500 m apart were used subsequently to select areas of hemodynamic data at 500 m intervals away from the electrode. This approach allowed us to investigate the spatial spread of positive and negative evoked hemodynamic responses in relation to the activated whisker barrel region. Statistical analysis. The data were statistically analyzed using two-tailed paired t tests, with Bonferroni-adjusted ␣ to correct for multiple comparisons. p levels for significance were set at 0.05 before Bonferroni's correction for multiple comparisons. For testing 16 s "train" stimulation trials, sorted by hemodynamic response magnitude, the sorted neural time series (top and bottom 33% of trials) were averaged for the period of stimulation (0 -16 s) and the t tests applied. Single-impulse stimulation and "nonstimulated" trials were sorted by the various neural candidates and the resulting trial-sorted hemodynamic time series (top and bottom 33% of trials) were averaged for the period 0 -5 s after stimulus presentation (equivalently, for nonstimulated trials, the same stimulus onset time and 0 -5 s poststimulation period were used even though a stimulus was not presented) and t tests applied.

Results
Neural responses were recorded from the activated whisker barrel region and surrounding regions that exhibited a negative hemodynamic response to sensory stimulation. The spatial and temporal dynamics of the stimulus-evoked hemodynamic responses are described below (see Stimulus-evoked positive and negative hemodynamic responses in the somatosensory cortex). Hemodynamic responses to 16 s trains of stimulation were sorted on an individual trial basis to evaluate their association with changes in MUA recorded from deep-layer regions (see Sorting negative hemodynamic responses from individual trials of 16 s train stimulation). Neural candidates were then investigated as potential markers of negative hemodynamic changes (see Sorting individual trials by neural activity from negative surround regions). Trial sorting by the long-latency gamma power was applied to the 2D hemodynamic changes evoked by the singleimpulse stimulation, with the spatial extent of the trial-sorted surround negative evaluated (see Spatial-temporal analysis of hemodynamic changes sorted by gamma power). To investigate the laminar profile of gamma power as a marker of negative hemodynamic magnitude, the long-latency gamma power from different cortical depths was used to sort the hemodynamics associated with each trial (see Laminar analysis of gamma-based power as a marker of the negative surround hemodynamic response). Finally, the long-latency gamma power was used to sort the 16 s train stimulation data and was compared with the trial-sorted single-impulse stimulation recordings (see Sorting individual 16 s stimulus duration (train-stimulation) trials by long-latency gamma activity from negative surround regions). In the interests of clarity and brevity, when reporting hemodynamic changes, Hbt is shown throughout the results section, as this aspect of the hemodynamic response is the composite of Hbr and HbO 2 and thus offers a robust measure of hemodynamic changes. Hbr and HbO 2 are included in overall representations of hemodynamic changes.

Stimulus-evoked positive and negative hemodynamic responses in the somatosensory cortex
The stimulus evoked increases in Hbt that were observed in the whisker barrel region (Fig. 1E, red regions to the lower right) were reliably associated with decreases in Hbt in surround regions (Fig. 1E, dark blue regions to the left of the whisker barrels). The consistent location of these positive and negative hemodynamic responses was used to guide the placement of the neural recording electrodes (Fig. 1G). Interestingly, clear positive hemodynamic responses were also observed in the whisker motor region following whisker pad stimulation (Fig. 1E, top left red regions showing increased Hbt).
Two stimulation paradigms (16 s stimulus train and singlepulse electrical stimulation) were used to evoke neural and ac-companying associated hemodynamic responses in the whisker barrel and surround regions of cerebral cortex. We found that the temporal dynamics of the hemodynamic responses associated with the different paradigms were dissimilar in both magnitude and profile. These are described respectively below.
Whisker stimulation: 16 s train Trains of 0.8 mA pad stimulation were presented for 16 s with an intertrial interval of 70 s ( Fig. 2A). The evoked neural activity recorded from the whisker barrel region produced an individual LFP depolarization to each of the 80 (5 ϫ 16) individual stimuli within the train (Fig. 2B), with a profile of gradually decreasing responses to successive stimuli. These electrophysiological responses were associated with a large-magnitude positive hemodynamic response recorded from the whisker barrel region (Fig.  2C). The hemodynamic changes showed an initial increase in Hbt, first peaking, before declining to a plateau phase, before returning to baseline ϳ10 s after stimulus cessation. These results confirm that trains of repeated stimuli, which are commonly used in block designs , produce large responses with increased signal-to-noise ratios. The neural activity evoked by whisker stimulation and recorded in regions surrounding the whisker barrel region was much reduced, with small changes in the LFP (Fig. 2B). The corresponding hemodynamic changes in the surrounding regions were also smaller in magnitude and had a negative profile, with a sustained decrease in Hbt for the duration of stimulus presentation, with an increase in Hbt following the cessation of stimulation (Fig. 2C, inset).

Whisker stimulation: single impulse
Single-impulse stimuli (300 s, 0.8 mA) were presented with a 25 s intertrial interval (Fig. 2D). The stimuli produced a robust and easily detected neural response (Fig. 2E) in whisker barrel cortex, while again in regions surrounding the whisker barrel's small-magnitude changes in LFP were observed in response to each stimulus. The hemodynamic responses to single stimuli ( Fig. 2F ) were much smaller than those evoked by the trains of electrical stimulation (Fig. 2C). The corresponding hemodynamic responses were also significantly smaller than those elicited by the 16 s stimulus train (Fig. 2,compare F,C). The positive hemodynamic responses evoked by single-pulse stimulation consisted of a small short-duration increase in Hbt, from which the following plateau phase was missing. Thus, after reaching an initial peak, the hemodynamic responses returned quickly to baseline with a poststimulus undershoot. Single-stimulus-evoked hemodynamic responses in regions surrounding the whisker barrels were also smaller but had a robust negative component (Fig. 2F,inset). This consisted of an initial short-duration decrease in Hbt, immediately followed by a poststimulus increase in Hbt.

Sorting negative hemodynamic responses from individual trials of 16 s train stimulation
Negative BOLD signals have been associated with decreases in both LFP and MUA (Shmuel et al., 2006). To investigate this phenomenon further, we ordered the stimulation trials according to the magnitude of the evoked negative hemodynamic response. We then applied the vector of sorted stimulation trials to concurrently recorded neural data. To compare trials that contained the largest and smallest negative responses, we averaged the top 33% and bottom 33% of trial blocks according to response magnitude (Fig. 3A). Note that the differences in hemodynamic responses between these two groups were due to spontaneous variation, as the stimulus conditions remained constant throughout.
Trial sorting of MUA by negative hemodynamic response Averaged, unsorted neural recordings from regions exhibiting negative hemodynamic responses showed increases in power at MUA frequencies in the superficial cortical layers (Fig. 3B, channels 2-4, equating to 150 -450 m below the cortical surface), but reductions in MUA power in the deeper cortical layers (Fig.  3B, channels 11-13, equating to 1050 -1350 m below the cortical surface), for the duration of stimulation. In the latter case there was a clear rebound increase in deeper-layer MUA power following stimulation offset. These results confirm our previously reported observations .
Sorting the electrophysiological data according to the magnitude of the concurrently recorded negative hemodynamic responses showed that trials with the smallest decreases in surround region Hbt had the greatest increase in shallow MUA power (Fig. 3C) and the smallest decrease indeep-layer MUA power (Fig. 3D). In contrast, trials that evoked the greatest surround region negative Hbt responses had the smallest increase in MUA power in superficial cortical layers, and the largest decrease in MUA power in the deep layers. Furthermore, trials with large negative hemodynamic responses tended to show larger rebound increases in deeplayer MUA power following stimulation offset. However, because the differences between the large and small negative hemodynamic response trials were clearest shortly after stimulus onset only, statistical analyses of MUA power covering the entire stimulation period revealed no reliable differences in either the superficial channel (paired t test: t ϭ 0.772, df ϭ 8, 3E). To uncover the neural drivers of the negative BOLD, we therefore turned our attention to the shorter-duration singleimpulse condition.

Sorting individual trials by neural activity from negative surround regions
Long-duration train stimulation paradigms ( Fig. 2A) are useful for evoking large-magnitude hemodynamic changes. However, the downside of this strategy is that long-latency neural responses and the complex interactions between individual evoked responses will be masked. Thus, our next step was to analyze the hemodynamic data after sorting individual stimulation trials by the neural activity evoked in negative surround regions by singleimpulse stimulation. The LFP and specific frequency bands of neural activity have been linked to negative BOLD (Shmuel et al., Boorman et al., 2010;Ramot et al., 2012). Thus, eight neural components were chosen for sorting of individual stimulation trials. The first being the stimulus-evoked LFP depolarization occurring in the whisker region and the remaining seven being based on frequency-band power of neural activity in the surround region. Our subsequent analyses sorted the hemodynamic data according to the magnitudes of each of the neural components.
Sorting neural recordings from individual stimulation trials is likely to be influenced by ongoing spontaneous neural activity. This phenomenon is exploited for functional mapping using resting-state recordings (Bruyns-Haylett et al., 2013). Therefore, in the present study we wanted to explore the contribution of spontaneous activity on individual trial sorting of neural recordings. To give us an estimate of the spontaneously occurring "natural variation" resulting from the sorting process, the hemodynamic data from nonstimulated "control" trials were sorted according to the magnitude of each of the eight neural components. Also, to confirm that our method of trial sorting did not impart inherent bias, control sorting was performed using randomized vectors of trials. When sorted in this manner, no significant differences were observed in the magnitudes of largest and smallest surround hemodynamic responses ( Fig. 4H; paired t test: t ϭ Ϫ0.0475, df ϭ 9, p ϭ 0.963, Bonferroni-corrected p ϫ 9 ϭ 8.669). This confirms that the trial-sorting procedure did not impart any inherent bias.

Single-impulse stimulation: sorting by LFPs
The LFP evoked by single-pulse whisker stimulation was easily distinguished with the main depolarization having a largemagnitude component lasting 10 -20 ms (Fig. 4A). Individual trials were therefore sorted into the largest 33% and smallest 33% according to the magnitude of the LFP occurring in the whisker barrel region (0 -20 ms after stimulus onset). These LFP-classified trials were then used to categorize corresponding hemodynamic data from the negative surround region (Fig. 4B). Stimulation trials with the smallest LFP depolarization had a slightly larger negative hemodynamic response, followed by a larger poststimulus positive component. However, comparing the hemodynamic responses (0 -5 s after stimulation) from negative surround regions of cortex that were associated with the largest and smallest whisker barrel LFP depolarizations (Fig. 4H ) revealed no reliable statistical differences (paired t test: t ϭ 0.564, df ϭ 9, p ϭ 0.587, Bonferroni-corrected p ϫ 9 ϭ 5.283). These results show the initial stimulus-evoked LFP depolarization in barrel cortex is unable to predict the magnitude of surround region negative hemodynamic responses.
Single-impulse stimulation: sorting by frequency-band power Next, individual trial recordings of neural activity from surround regions that exhibited negative hemodynamic responses were subjected to frequency-power analysis. Neural activity was first placed in 100 ms temporal bins before applying Fourier transforms (Fig. 4C,F ), thus allowing observation of the temporal dynamics of the neural activity within the selected frequency bands. Note that the initial increase in power (0 -100 ms) observed across all channels in both MUA (300 -3000 Hz) and gamma (30 -80 Hz) frequency bands is unreliable as this period contains the initial LFP depolarization, which has a nonsinusoidal shape, ensuring that the frequency spectra of this signal are very broadband. This produces aliasing and leakage across frequencies when basic Fourier analysis techniques are applied. The 0 -100 ms period following stimulus onset was therefore excluded from trial sorting using frequency-power analysis. For trial sorting by the various neural candidates, activity was not placed into 100 ms bins, but was instead taken in continuous form from the selected time ranges and subject to Fourier analysis to improve both the frequency resolution and quality of the frequency-power analysis. Decreases in MUA power have previously been observed in the deep cortical layers of regions exhibiting sustained negative hemodynamics evoked by trains of sensory stimuli . Therefore, the hemodynamic responses in surround regions were grouped according to trials with the 33% largest and 33% smallest changes in MUA power recorded from surround region deep cortical layers (1000 -1300 m below the cortical surface; Fig. 4E). The analysis was restricted to the period 100 -1000 ms after stimulus presentation to match the initial singleimpulse neural response (Fig. 2E). Comparison of surround hemodynamic data sorted on the basis of local deep-layer MUA power revealed modest differences (Fig. 4E), but failed to achieve statistical significance (paired t test: t ϭ 2.262, df ϭ 9, p ϭ 0.050, Bonferroni-corrected p ϫ 9 ϭ 0.450). This finding shows that deep-layer MUA power is another unreliable correlate of the main negative component of the surround hemodynamic changes evoked by single-impulse stimuli.
Increases in MUA power (0 -300 ms from stimulus onset) were observed in the upper cortical layers (400 -900 m below the cortical surface) of the negative surround region following the single-impulse stimuli (Fig. 4C). This region is known to have major input from the thalamus (Armstrong-James et al., 1992). The individual-trial hemodynamic data from the negative surround regions were therefore sorted according to the 33% largest and 33% smallest changes in MUA power recorded from the superficial layers (250 -750 m below the cortical surface; Fig.  4C). The "long-latency" time range, from 300 to 2000 ms after stimulus presentation, was chosen to include a period of decreased MUA power observed following stimulus presentation in the whisker barrel region (data not shown), while avoiding the large increases in MUA power immediately following the stimulus (Fig. 4C, black box). The hemodynamic responses in trials with high and low long-latency MUA power in the superficial layers were significantly different (Fig. 4D). A small transient increase in Hbt was observed for 33% of trials with higher than average MUA power, while a larger decrease in Hbt was observed for the 33% of trials having the least MUA power (paired t test: t ϭ 4.294, df ϭ 9, p ϭ 0.0020, Bonferroni-corrected p ϫ 9 ϭ 0.018). These results show that long-latency (300 -2000 ms poststimulus) MUA power in the upper cortical layers is predictive of surround region negative hemodynamic responses.
To discover whether there is a better neural predictor of negative hemodynamic changes, we turned our attention to the lower-frequency LFP oscillations: delta-band, theta-band, alphaband, beta-band, and, in particular, gamma-band power, which has been linked to BOLD in several previous studies (Lachaux et al., 2007;Nir et al., 2007;Goense and Logothetis, 2008;Ramot et al., 2012). In the current trial-averaged neural recordings, we observed an increase in gamma power that lasted for ϳ600 ms after stimulus presentation and was concentrated in the upper cortical layers of the surround region (Fig. 4F ). The singleimpulse stimulation trials were therefore sorted by gamma power between 250 and 750 m below the cortical surface, 300 -2000 ms after stimulation (Fig. 4F, black box). Sorting by long-latency gamma power produced large differences between the time series of the surround region hemodynamics (Fig. 4G). The 33% of trials having the most gamma power were associated with an increase in Hbt, while the 33% of trials having the least longlatency gamma power had a large decrease in Hbt following stimulation ( Fig. 4H; paired t test: t ϭ 3.675, df ϭ 9, p ϭ 0.0051, Bonferroni-corrected p ϫ 9 ϭ 0.0459). Shown for comparison are the results for the other common EEG bands (Fig. 4H ). Sorting individual hemodynamic trials by the delta and theta frequency bands did not result in significant differences between the hemodynamics responses in the 33% of trials having the most and least power in these bands (delta: paired t test, t ϭ 0.176, df ϭ 9, p ϭ 0.862, Bonferroni-corrected p ϫ 9 ϭ 7.761; theta: paired t test, t ϭ Ϫ0.161, df ϭ 9, p ϭ 0.875, Bonferroni-corrected p ϫ 9 ϭ 7.879). Observation of the 33% of trials, sorted by the most and   least alpha and beta power, suggested they hemodynamic responses differed in magnitude, but statistical analysis did not reach significance (alpha: paired t test, t ϭ 2.347, df ϭ 9, p ϭ 0.0435, Bonferroni-corrected p ϫ 9 ϭ 0.392; beta: paired t test, t ϭ 3.061, df ϭ 9, p ϭ 0.0136, Bonferroni-corrected p ϫ 9 ϭ 0.122). However, the largest difference in the magnitude of the surround region hemodynamics was observed when sorting trials by the long-latency gamma band. Thus, these results show that long-latency gamma offers a robust colocalized prediction of negative hemodynamic changes, outperforming the other chosen neural sorting candidates (Fig. 4H ).
To confirm that the responses illustrated in Figure 4 were attributable to single-pulse stimulation of the whisker pad, data were collected using the same paradigm (100 trials, with 25 s ISI) in the absence of any electrical stimulus. Individual "nonstimulation" trials were sorted in an identical manner to that used to sort the single-impulse stimulation data. Each of the nine neural candidates (initial LFP depolarization, deep MUA, shallow MUA, and the classic EEG frequency bands) were again used to sort the nonstimulated data. The trial-sorted hemodynamics (top and bottom 33% of trials) were compared by taking the averaged hemodynamic response for the period 0 -5 s after the time point where a stimulus would have been presented in the stimulation data. No reliable differences were observed in the hemodynamic responses associated with the top and bottom 33% of trials when the data were sorted by any of the neural candidates (paired t test: initial LFP depolarization, t ϭ Ϫ5.096, df ϭ 4, p ϭ 0.0070, Bonferroni-corrected p ϫ 8 ϭ 0.0560; MUA deep, t ϭ 0.635, df ϭ 4, p ϭ 0.560, Bonferroni-corrected p ϫ 8 ϭ 4.478; long-latency, MUA shallow, t ϭ 1.3894, df ϭ 4, p ϭ 0.237, Bonferronicorrected p ϫ 8 ϭ 1.897; gamma, t ϭ 1.297, df ϭ 4, p ϭ 0.265, Bonferroni-corrected p ϫ 8 ϭ 2.116; delta, t ϭ Ϫ0.973, df ϭ 4, p ϭ 0.386, Bonferroni-corrected p ϫ 8 ϭ 3.087; theta, t ϭ 0.137, df ϭ 4, p ϭ 0.898, Bonferroni-corrected p ϫ 8 ϭ 7.181; alpha, t ϭ Ϫ0.292, df ϭ 4, p ϭ 0.785, Bonferroni-corrected p ϫ 8 ϭ 6.276; beta, t ϭ 2.754, df ϭ 4, p ϭ 0.0511, Bonferroni-corrected p ϫ 8 ϭ 0.409). Thus, these results show that spontaneous activity was unlikely to have had any significant influence on the hemody-namic responses when trials were sorted by stimulation-evoked neural activity.
The change in HbO 2 sorted by long-latency gamma showed even greater differences between the hemodynamic time series (Fig. 5 B, C). However, it is the time series of Hbr that shows the most noteworthy change, being the most closely linked to BOLD. The greatest increase in Hbr can be observed for stimulation trials having the least gamma power (Fig. 5B) and it is this signal that would generate the largest negative BOLD response. Conversely, the trials having the most gamma power showed no substantial change in Hbr immediately following stimulation, and are therefore unlikely to generate a negative BOLD signal (Fig. 5C).

Spatial-temporal analysis of hemodynamic changes sorted by gamma power
Having established long-latency gamma power as the best predictor of associated hemodynamic response, we next investigated the spatial-temporal nature of the evoked hemodynamic changes in the surround region. A representative example of 2D-OIS recordings of hemodynamic changes sorted by long-latency gamma power is illustrated in Figure 5D-G. The spatially averaged hemodynamic changes from all trials showed an increase in Hbt overlying the whisker barrels observable 2-3 s following stimulus presentation ( Fig. 5E; Movie 1), with a surrounding small decrease in Hbt. Approximately 6 s after stimulus presentation, the pattern was reversed with a decrease in Hbt in the whisker area accompanied by an increase in the surround region.
The trial-averaged hemodynamic changes with the least longlatency gamma power showed the greatest decrease in Hbt 2-3 s after stimulus presentation in the region surrounding the whisker barrels ( Fig. 5F; Movie 2), and the greatest subsequent increase ϳ4 s later. Thus, the stimulus sets up opposed but decaying oscillatory hemodynamic responses in the barrel and surround regions with the peaks of responses in the surround region delayed by 0.5-1.0 s ( Fig. 5F; Movie 2). The hemodynamic changes from the trials having the most long-latency gamma power showed little change in Hbt in surround regions ( Fig. 5G; Movie 3). These results show that long-latency gamma power exposes differences in the spatial-temporal nature of the hemodynamic response in the negative surround region.
To estimate the spread of hemodynamic responses from the whisker barrel region across all subjects, a concentric ring analysis was developed (Fig. 5H-K ). Centered on the whisker barrel cortex of each subject, 10 concentric rings (increasing by ϳ500 m) were overlaid onto the 2D optical imaging data so that the outer rings encompassed the area of maximum surrounding Hbt decrease (Fig. 5H ). For each subject, pixels inside each ring were averaged to produce a time series of Hbt change. Concentric ring plots were generated for group-averaged data and the datasets sorted by the long-latency gamma-band neural activity (Fig. 5I-K ). Stimulation trials having the least long-latency gamma power in the surround region showed a larger-magnitude spatially defined negative region, ϳ2-3 mm from the whisker barrel electrode (Fig. 5J ). Trials having the most long-latency gamma power in the surround regions showed a minimal decrease in Hbt (Fig.  5K ). A stimulus-evoked oscillation was observed in both the average of all trials and in trial-sorted datasets (Fig. 5I-K ). Responses in the surrounding region were phase lagged with respect to those occurring in the whisker region. These results illustrate the complex spatial nature of the surround hemodynamic H, Comparison of trial sorting by each neural candidate. Each bar represents the change in Hbt in the surrounding regions averaged between 0 and 5 s after stimulus presentation, with the top and bottom 33% of trials selected dependent on the various neural markers or a randomized vector of trials. The green dotted line represents the surround region Hbt response averaged between 0 and 5 s after stimulus presentation and averaged across all stimulation trials. The line offers a reference to show the difference from the mean of each of the Hbt responses that arise after sorting by each neural metric. The overlaid blue box represents the hemodynamic trials sorted by each frequency range, taking the neural data for the long latency period after stimulus presentation (300 -2000 ms) and from channels 3-8 of the surround electrode. Paired t tests used to assess reliability of each sorting approach. *p Ͻ0.05 (Bonferroni corrected for multiple comparisons ␣ ϭ 0.05/9, p level Ͻ0.0056). Error bars represent SEM across subjects. Responses averaged across subjects (n ϭ 10). changes, which appear to interact with positive responses in the activated whisker barrel region.

Laminar analysis of gamma-band power as a marker of the negative surround hemodynamic response
To make full use of the depth resolution of the 16-channel linear array electrodes, trial sorting was performed on a channel-bychannel basis. The averaged hemodynamic responses from the top and bottom 33% of trials sorted by gamma power were taken for each channel. The overall mean from all trials was subtracted to align the effects of sorting by long-latency gamma about the zero axes. The change in Hbt between 0 and 5 s after stimulus presentation was calculated for each channel for both stimulus and spontaneous (nonstimulation control) trials (Fig. 6). The stimulus-evoked trials having least gamma power showed the greatest consistency across subjects, with the channels centered at 400 m below cortical surface producing the greatest negative hemodynamic response. Stimulation and spontaneous control trials having most gamma-band activity showed highly variable increases in surround Hbt across all 16 channels (Fig. 6). To confirm the laminar specificity of the stimulus-evoked gammaband trial sorting, the two original cortical depth ranges used for sorting (Fig. 4), 250 -750 m and 1050 -1350 m, were compared (Fig. 6, black and pink boxes). A one-way ANOVA (F (1,18) ϭ 14.05; p ϭ 0.0015) confirmed that the decrease in long-latency gamma power in the upper cortical layers, peaking at ϳ400 m below the cortical surface, is best able to "predict" stimulusevoked negative hemodynamic changes.

Sorting individual 16 s stimulus duration ("train-stimulation") trials by long-latency gamma activity from negative surround regions
The surround region hemodynamics generated by 16 s train stimulation data were used to sort the concurrently recorded MUA (see Sorting negative hemodynamic responses from individual trials of 16s train stimulation). However, there were not any reliable differences in MUA in the deep or shallow layers. Because the sorting of single-impulse stimulation trials by the long-latency gamma-band power produced robust differences in the surround region hemodynamic changes, the same analytical  Figure 5. Spatial-temporal representation of hemodynamic changes evoked by single-impulse stimulation, sorted by gamma power (30 -80 Hz) from the upper cortical layers (channels 3-8) of the surround electrode (300 -2000 ms after stimulus presentation). A, Hemodynamic response in surround region, evoked by single-impulse stimulation of the whisker pad, average of all trials. B, Hemodynamic response in surround region, average of trials (33%) having the least gamma power in the surround region, in the upper cortical layers (channels 3-8). C, Hemodynamic response in surround region, average of trials (33%) having the most gamma power in the surround region, in the upper cortical layers. D, In vivo grayscale image of thinned cranial window from a single representative subject, with regions of interest around the two recording electrodes used for time-series generation. E-G, Averaged Hbt changes, evoked by single-impulse stimulation from a single representative subject, all trials (100%), least gamma (33%) and most gamma (33%) respectively, and gamma power recorded from surround electrode. Hbt changes in the whisker region response (black) and surround region response (green). Montage of images shows spatial maps of micromolar changes in Hbt. Each is averaged over a second and represents a corresponding time point in the stimulation trial. Increases in Hbt are toward the red color range, with blue representing negative changes. H, In vivo grayscale camera image of thinned cranial window from a single representative subject, with a series of concentric rings centered on the whisker barrel electrode. These were used to select pixels for selection of data as a function of distance away from the electrode. I-K, Hbt changes, with data from all trials, least gamma power (33%) and most gamma power (33%) respectively. The y-axis represents distance away from the tip of the whisker electrode, with each row in the image generated from pixels within each concentric ring. Stimulation represented with arrows. Error bars represent SEM across subjects. Responses averaged across subjects (n ϭ 10), except for single subject (D-H).
approach was again applied to the 16 s train stimulation data. The gamma-band (30 -80 Hz) frequency power analysis was applied to the 16 s train stimulation recordings, with data placed into 100 ms bins (Fig. 7A). An increase in gamma power was observed across all cortical layers for the duration of the 16 s train stimulation, with a localized increase during each individual stimulus pulse of the 16 s stimulation train. The individual trials were sorted by again taking the gamma-band power from electrode channels 3-8, equating to 250 -750 m below the cortical surface. However, two time ranges were used for sorting. The first used the original range previously used for sorting the singleimpulse stimulation recordings (300 -2000 ms; Fig. 7A, black box). The second range was matched to the longer components of the 16 s train stimulation (300 -18,000 ms; Fig. 7A, pink box). Statistical analysis was performed in the same manner as for the single-impulse stimuli, using two-tailed paired t tests, with Bonferroni-adjusted ␣ to correct for multiple comparisons.
Sorting the surround region hemodynamic changes, evoked by the 16 s train stimulation, by the long-latency gamma-band power (300 -2000 ms) did not produce clear observable differ-ences between the 33% of trials having the most gamma and the 33% of trials having the least gamma (Fig. 7B). To assess whether sorting by the long-latency gamma could predict the magnitude of the initial surround hemodynamic response, the trial-sorted and averaged responses were averaged between 0 and 5 s after stimulus onset (Fig. 7D) and were subjected to statistical analysis. Although the trials having the least long-latency gamma had the most negative average Hbt changes between 0 and 5 s after stimulation, there were no significant difference between the 33% of trials having the most and the 33% of trials having the least gamma (paired t test, t ϭ 1.0828, df ϭ 5, p ϭ 0.3283, Bonferronicorrected two comparisons p ϫ 2 ϭ 0.6566). The sorting of the surround region hemodynamic changes by the extended period of gamma (300 -18,000 ms after stimulation onset) was also unable to produce any clearly observable differences between the 33% of trials having the most and the 33% of trials having least gamma power (Fig. 7C). The Hbt time series (33% of trials having the most gamma and 33% of trials having the least gamma) were averaged between 0 and 21 s after stimulus onset to assess the ability of sorting by gamma over the whole period of stimulation (Fig. 7D). Again the trials having the least gamma had the mostnegative average Hbt between 0 and 21 s after stimulation. However, the difference between the most and least gamma trials was not significant (paired t test, t ϭ 1.2613, df ϭ 5, p ϭ 0.2628, Bonferroni-corrected two comparisons p ϫ 2 ϭ 0.5256). This shows the long-latency gamma marker is unable to reliably predict the negative hemodynamic response generated by longduration trains of stimuli.

Discussion
Negative hemodynamic responses evoked by single-impulse stimulation of the whisker pad are closely associated with longlatency reductions in gamma-band (30 -80 Hz) neural activity in somatosensory cortex adjacent to activated whisker barrel cortex. The magnitude of the negative hemodynamic response was predicted best by decreases in gamma activity in middle-superficial cortical layers (II-IV). Spatiotemporal analysis of hemodynamic responses indicated that single-impulse stimuli initiated a decaying series of spatially distributed positive and negative hemodynamic oscillations in which the negative response in surrounding regions had a partially delayed antiphase compared with the positive response in the barrel region. These results suggest that neg-Movie 1. Micromolar change in Hbt for a single representative animal, average responses from all stimulation trials. Lower right, Grayscale CCD camera image of thinned cranial window overlying whisker barrel field, with overlaid whisker (black) and surround regions of interest (green). Lower left, Hbt change generated from regions of interest.
Movie 2. Change in Hbt for a single representative animal, average response from 33% of stimulation trials having the least long-latency (300 -2000 ms after stimulus presentation) gamma-band activity from surround electrode channels 3-8.
Movie 3. Change in Hbt for a single representative animal, average response from 33% of stimulation trials having the most long-latency (300 -2000 ms after stimulus presentation) gamma-band activity from surround electrode channels 3-8.

Micromolar change in Hbt
Evoked ( ative BOLD signals reflect genuine changes in neural activity, rather than a redistribution of blood ("vascular steal") into the neighboring activated region (Harel et al., 2002). Consistent with our previous work , the negative BOLD responses evoked by long-duration trains of stimulating pulses were associated with decreases in MUA in deep cortical layers adjacent to the activated whisker barrel field. Importantly, single-impulse sensory stimuli failed to elicit the same phenomenon. In a similar manner, the long-latency gamma-band power was able to predict negative hemodynamic changes evoked by single-impulse stimuli, but not for longerduration 16 s train stimuli. These differences could arise from the interactions of long-latency neural activity between successive pulses of the stimulus train. For example, inspection of the gamma-power changes for the single-impulse stimuli (Fig. 4F ) shows an increase in gamma across all layers, decaying 300 ms after stimulus presentation, while around layer IV the increase in gamma lasts over 600 ms after stimulus presentation. These gamma increases are occurring over periods longer than the spacing between successive trains of stimulus pulses (5 Hz corresponds to 200 ms). This disparity is most likely explained by the compounding effects of repeated stimuli that induce sensory adaptation and nonlinearities in the stimulus-response relationship. These could conceal or reveal different processes influencing the responses to successive stimuli presented at different time scales. This highlights the importance of using single stimuli to investigate the basic mechanisms underpinning neurovascular coupling.
Surround negative hemodynamic responses evoked by somatosensory stimulation have been suggested by some to be a "capricious phenomenon" (Lauritzen et al., 2012). However, in this study we have shown that sustained, regionally specific negative hemodynamic responses were elicited reliably both within and across subjects. Consistent negative hemodynamic responses were observed over the forepaw, hindpaw, and trunk cortical regions that are medially adjacent to the activated whisker barrel field. Similar results in human visual cortex using fMRI have been reported (Smith et al., 2010). We propose that the distinctive spatial localization of the negative hemodynamic response could reflect a mechanism for suppressing somatosensory representations with less stimulus relevance than the activated whisker barrel field. This suggestion is in keeping with a recent study investigating negative BOLD in humans that proposed that changes in gamma-band activity represented a widespread mechanism for modulating cortical processing (Ramot et al., 2012).
Gamma oscillations in cerebral cortex are thought to be generated by the synchronous firing of multiple inhibitory GABAergic interneurons (Traub et al., 1998;Cardin et al., 2009). They are believed to reflect higher cognitive processing, including memory formation (Jutras and Buffalo, 2010) and sensory perception (Lachaux et al., 2005). Despite much research on the neurophysiological basis of gamma oscillations, little is known about their functional role in information processing (Bartos et al., 2007) or indeed their link to the neurovascular coupling (Riera and Sumiyoshi, 2010).
Converging evidence suggests that increases in gamma-band activity are closely associated with positive hemodynamic responses in health (Niessing et al., 2005) and disease (Harris et al., 2014). Moreover, they seem to be strongly correlated with both resting-state and stimulus-evoked BOLD fMRI signals (Lachaux et al., 2007;Nir et al., 2007;Goense and Logothetis, 2008). While these studies are consistent with the view that BOLD signaling should be positively related to higher-frequency EEG activity (Kilner et al., 2005), it has been unclear whether decreased gamma-band activity is coupled with negative BOLD responses (Lachaux et al., 2007;Yuan et al., 2011). Concentrations of inhibitory GABA in human anterior cingulate cortex predict negative BOLD responses (Northoff et al., 2007), and resting state GABA concentrations correlate inversely with positive BOLD fMRI responses (Muthukumaraswamy et al., 2009). Multiple pathways have been proposed by which increased gamma firing can produce local vasodilation leading to positive hemodynamic changes (Riera and Sumiyoshi, 2010). Reductions in gamma-band activity could therefore produce either active or passive vasoconstriction by modulating similar pathways in various ways, such as by reducing nitric oxide release or increasing arachidonic acid (Attwell et al., 2010). An additional explanation for the association between the changes in gamma power and hemodynamics could be related to the metabolic requirements of gamma oscillation. The need to maintain ionic homeostasis during gamma oscillations has been shown to be metabolically demanding (Kann et al., 2011;Lord et al., 2013). A reduction in gamma power could therefore be associated with reduced metabolic demand, leading to a local vasoconstriction and a subsequent negative hemodynamic response.
Inhibitory interneurons account for a small proportion (ϳ20 -30%) of neurons in rodent neocortex. Approximately half of these are basket cells that project to pyramidal output neurons and other cortical interneurons. These inhibitory interneurons have lateral projections that span cortical columns (Markram et al., 2004). Cortical layers II-IV contain the greatest density of inhibitory interneurons (Markram et al., 2004;Meyer et al., 2011). These layers exhibit strong inhibitory drive, so that the columnar spread of excitatory activity from layer IV to layer II/III engages more inhibitory than excitatory neurons (Mateo et al., 2011). Preferential connectivity between excitatory neurons in all cortical layers and layer II/III inhibitory interneurons has also been demonstrated (Dantzker and Callaway, 2000). Afferent inhibitory projections from the thalamus that terminate in cortical layers II/III (Timofeev andChauvette, 2011) andIV (Porter et al., 2001) can provide external inhibitory input to supplement intrinsic cortical inhibition. These thalamic relay neurons operate either in burst or tonic modes (Jahnsen and Llinás, 1984;Sherman, 2001) and could be an important modulating source of cortical oscillatory activity (for review, see Jones 2009). Our finding that a long-latency gamma-band neural marker of negative hemodynamic responses was particularly prevalent in cortical layers II-IV is consistent with the above research and provides new insights into the important role played by these laminae in cortical function and resultant hemodynamics.
The decrease in gamma-band activity 300 -2000 ms after stimulus presentation was the neural parameter that best predicted the magnitude of the evoked negative hemodynamic change. This delayed gamma response appears too protracted to be explained by local operations within cortical circuits, e.g., direct modulation by stimulus-evoked activity in adjacent barrel cortex. Rather, it seems more likely that the long-latency decrease in gamma activity could be mediated through cortical interactions with subcortical circuits. One potentially relevant structure would be the thalamic reticular nucleus (TRN), which is considered to be an important driver of inhibition within the whisker barrel field (Temereanca and Simons, 2004). The TRN possesses rich interconnectivity between constituent nuclei and has topographically mapped connections with cortex that enable both discrete and diffuse inhibitory control (Guillery et al., 1998). Additionally, the TRN has been implicated in the generation of gamma oscillations in auditory cortex (Macdonald et al., 1998) and, in slice, has been shown to respond differentially to stimulation at gamma and nongamma frequencies (Mistry et al., 2008). Moreover, TRN neurons both in anesthetized and awake rats show slow recovery from excitation (Yu et al., 2009). Thus, auditory stimuli were found to produce diminished neural responses in the TRN when presented at intervals between 200 and 1000 ms and were found to have fully recovered when stimuli were presented at intervals over 3000 ms. This recovery period is within the 2000 ms poststimulus window of the long-latency decrease gamma associated with the negative hemodynamic response found in the present study.
An alternative mechanism to account for the delayed gamma response could involve the re-entrant looped architecture connecting somatosensory cortex with the basal ganglia. Excitatory input to somatosensory cortex from the thalamus (including the ventromedial nucleus) is kept under inhibitory control by the basal ganglia output nuclei (substantia nigra pars reticulate and the entopeduncular nucleus; Chevalier and Deniau, 1990). Increased inhibitory control of thalamic projections to surrounding nonbarrel somatosensory cortex, and selective disinhibition of loops returning to the activated barrel cortex could produce the appropriate cortical modulations at the required latencies following stimulus onset (Redgrave et al., 1999). In both humans and animals, the basal ganglia have been linked to selective attention (Kinomura et al., 1996). Dysfunction within this system has been associated with decreases in attentional performance (De-Long and Wichmann, 2007;Smith et al., 2011) and sensorimotor neglect (Marshall et al., 1974). Consequently it may be significant that gamma oscillations are attention dependent and are modulated in neuropsychiatric disorders associated with basal ganglia dysfunction, including schizophrenia and Alzheimer's disease (for review, see Herrmann and Demiralp, 2005). The anatomical and functional properties of the thalamic reticular nucleus and basal ganglia suggest these subcortical circuits may be promising candidates for mediating the long-latency surround inhibition in somatosensory cortex observed in the present study.
The use of anesthesia in the present study was necessary to enable the prolonged acquisition of high-resolution multimodal data using invasive techniques. Although a range of anesthetic regimes are used in neurovascular research, urethane provides for a highly stable, long-duration state of anesthesia that has repeatedly facilitated advances in our understanding of neurovascular coupling (Devor et al., 2007;Boorman et al., 2010). Notwithstanding this, urethane, like all anesthetic regimes, has the potential to disrupt normal neuronal, neurovascular, and hemodynamic response mechanisms (for a full discussion, see recent reviews by Masamoto and Kanno, 2012;Martin, 2014). However, the negative-surround hemodynamic response used here as a model for our neurophysiological investigation of negative BOLD signaling has been robustly observed under other anesthetic conditions (Devor et al., 2007) as well as in awake rats (Martin et al., 2013). Furthermore, urethane anesthesia has been shown to provide long-term stability and balanced actions on multiple neurotransmitter receptors (Masamoto and Kanno, 2012), while its neurovascular coupling actions appear to be mediated via changes in cortical excitability and slight alterations in the timing of the hemodynamic response function (Martin et al., 2006), neither of which seem likely to be able to account for the negative BOLD signaling relationships we observe here. Verification of the key findings reported here in other laboratories using alternative anesthetic regimes will, however, be an important next step.