Evidence that the negative BOLD response is neuronal in origin: A simultaneous EEG–BOLD–CBF study in humans

Unambiguous interpretation of changes in the BOLD signal is challenging because of the complex neurovascular coupling that translates changes in neuronal activity into the subsequent haemodynamic response. In particular, the neurophysiological origin of the negative BOLD response (NBR) remains incompletely understood. Here, we simultaneously recorded BOLD, EEG and cerebral blood flow (CBF) responses to 10 s blocks of unilateral median nerve stimulation (MNS) in order to interrogate the NBR. Both negative BOLD and negative CBF responses to MNS were observed in the same region of the ipsilateral primary sensorimotor cortex (S1/M1) and calculations showed that MNS induced a decrease in the cerebral metabolic rate of oxygen consumption (CMRO2) in this NBR region. The ∆CMRO2/∆CBF coupling ratio (n) was found to be significantly larger in this ipsilateral S1/M1 region (n=0.91±0.04, M=10.45%) than in the contralateral S1/M1 (n=0.65±0.03, M=10.45%) region that exhibited a positive BOLD response (PBR) and positive CBF response, and a consequent increase in CMRO2 during MNS. The fMRI response amplitude in ipsilateral S1/M1 was negatively correlated with both the power of the 8-13 Hz EEG mu oscillation and somatosensory evoked potential amplitude. Blocks in which the largest magnitude of negative BOLD and CBF responses occurred therefore showed greatest mu power, an electrophysiological index of cortical inhibition, and largest somatosensory evoked potentials. Taken together, our results suggest that a neuronal mechanism underlies the NBR, but that the NBR may originate from a different neurovascular coupling mechanism to the PBR, suggesting that caution should be taken in assuming the NBR simply represents the neurophysiological inverse of the PBR.


Introduction
Blood oxygenation-level dependent (BOLD) functional magnetic resonance imaging (fMRI) is widely used for non-invasive measurement of the spatial location and intensity of human brain activity. An increase in neuronal activity in response to stimulation results in a local increase in blood oxygenation and a corresponding increase in BOLD signal (i.e. a positive BOLD response (PBR)) relative to the prestimulus baseline period (Ogawa et al., 1990). The PBR amplitude has been shown to correlate most strongly with increases in local field potential (LFP) measurements of neuronal activity (Logothetis et al., 2001;Magri et al., 2012;Viswanathan and Freeman, 2007). A decrease in BOLD signal below pre-stimulus baseline levels is termed a negative BOLD response (NBR). A NBR has been observed to occur spatially adjacent to the PBR in visual cortex (Bressler et al., 2007;Pasley et al., 2007;Shmuel et al., 2002;Smith et al., 2004;Tootell et al., 1998;Wade and Rowland, 2010). In addition, a NBR has also been recorded in the motor (Allison et al., 2000;Hamzei et al., 2002;Newton et al., 2005;Stefanovic et al., 2004) and somatosensory (Hlushchuk and Hari, 2006;Kastrup et al., 2008;Klingner et al., 2010Klingner et al., , 2011a cortices, ipsilateral to the presented stimulus, with a concurrent PBR observed in the contralateral hemisphere. The ability to measure and create functional maps of neuronal deactivation using the NBR could provide important insights into the functional role of inhibition in brain processes. However, caution is required when interpreting decreases in the BOLD signal, as BOLD contrast arises from a complex interaction between changes in cerebral blood flow (CBF), cerebral blood volume (CBV) and cerebral metabolic rate of oxygen consumption (CMRO 2 ). Whilst the exact physiological origin of the NBR remains poorly understood, several generating mechanisms have been proposed: 2. A reduction in CBF due to a decrease in neuronal activity with a lesser reduction in CMRO 2 , resulting in an increase in the oxygen extraction fraction (OEF) (Pasley et al., 2007;Shmuel et al., 2002); 3. An increase in neuronal activity and CMRO 2 without a compensatory increase in CBF, again leading to an increase in the OEF (Schridde et al., 2008); 4. A CBF decrease in upper cortical layers accompanied by a CBV increase in middle cortical layers (as demonstrated in the visual cortex by (Goense et al., 2012)).
These conflicting findings make it difficult to interpret the NBR unambiguously without additional physiological measurements. Such measurements can be provided by multi-modal imaging.
Recent studies have shown that a haemodynamic steal cannot fully explain the NBR (Chen et al., 2005;Kastrup et al., 2008;Smith et al., 2004) with conservative estimates suggesting that at least 60% of the NBR is due to a neural component (Pasley et al., 2007;Shmuel et al., 2006). Invasive recordings demonstrated a decrease in both broadband LFP and spiking activity in NBR regions of primate visual cortex (Shmuel et al., 2006) and of somatosensory cortex in rats (Boorman et al., 2010). In humans, a decrease in CBF and CMRO 2 has been observed to occur in NBR regions of both the motor (Stefanovic et al., 2004) and visual (Shmuel et al., 2006) cortex. The average magnitude of the NBR has been shown to increase with increasing stimulus duration and intensity (Klingner et al., 2010;Shmuel et al., 2002), analogous to the behaviour of the PBR, suggesting that the PBR and NBR share a common link to the neuronal activity that is induced by stimulation. These studies suggest that a reduction in CBF, combined with a lesser reduction in CMRO 2 (mechanism 2), is the most likely physiological cause of the NBR. Increasing evidence suggests that the NBR reflects a functionally relevant measure of neuronal deactivation (Bressler et al., 2007;Ferbert et al., 1992;Kastrup et al., 2008;Klingner et al., 2010;Liepert et al., 2001;Schafer et al., 2012). However detailed clarification of the relationship between NBR, CBF and measures of the underlying neuronal activity in humans is still needed.
The simultaneous recording of electroencephalography (EEG) during BOLD fMRI allows the integration of brain activity measurements from neuronal and haemodynamic sources, and has been used to investigate correlations between the amplitude of neuronal responses (e.g. evoked potentials, alpha, beta or gamma oscillatory activity) and the PBR in sensory and cognitive studies (Debener et al., 2005;Eichele et al., 2005;Mayhew et al., 2010a;Mobascher et al., 2009;Mulert et al., 2010;Novitskiy et al., 2011;Ritter et al., 2009;Sadaghiani et al., 2010;Warbrick et al., 2009;Yuan et al., 2010). These studies have demonstrated that simultaneous EEG-fMRI can provide greater specificity regarding the spatial arrangement (Goldman et al., 2009;Novitskiy et al., 2011) or temporal sequence (Eichele et al., 2005) of responses in brain areas, compared to that revealed by a standard analysis of data from a single neuroimaging modality. The implementation of simultaneous arterial spin labelling (ASL) in conjunction with EEG-BOLD recordings provides additional measurements of changes in CBF, a major physiological component of the BOLD signal. ASL provides increased sensitivity to haemodynamic effects in arterial-capillary vessel network compared with BOLD contrast which is more venous in nature, thus its origin is more localised to regions of neuronal activity (Duong et al., 2001) and may therefore exhibit greater sensitivity to natural variations in this activity between trials. Simultaneous CBF-BOLD recordings also provide the possibility of evaluating changes in CMRO 2 within NBR regions, under the assumption that a similar coupling between BOLD, CBF, CBV and CMRO 2 occurs in regions of NBR as is observed in PBR areas (Shmuel et al., 2002;Stefanovic et al., 2004). Therefore EEG-BOLD-ASL measurements provide a powerful non-invasive imaging strategy for probing neurovascular coupling (e.g. (Mullinger et al., 2013b)) and as such may allow differentiation between the different potential origins of the NBR in humans.
The aim of the current study is to investigate the relationship between simultaneously recorded NBR, CBF and EEG responses to consistent intensity median nerve stimulation (MNS) applied to the right wrist. Employing lateralised MNS allows the NBR to be studied in ipsilateral sensorimotor cortex (S1/M1) whilst minimising the potential confound of blood steal by the PBR in contralateral S1/M1, due to the distinct vascular territories of the two cortical hemispheres (Tatu et al., 1998(Tatu et al., , 2012. Furthermore, unilateral MNS allows the NBR to be interrogated in terms of two types of EEG responses: evoked potentials which are transient changes in the ongoing EEG signal, phase-locked to stimulation (Hari et al., 1993;Kakigi, 1994); and non phase-locked event-related synchronisation/desynchronisation (ERS/ERD) of alpha (8-13 Hz) and beta (13-30 Hz) frequency oscillatory activity (Hirata et al., 2002;Neuper et al., 2006;Stevenson et al., 2012). Here, we use the general linear model (GLM) to investigate trial-by-trial correlations between natural variations in these EEG signals and both BOLD and CBF responses.
By combining estimates of oxygen metabolism (CMRO 2 ) with CBF and EEG measurements we aim to shed new light on the origin of the negative BOLD signal in humans. We hypothesise that the amplitude of BOLD responses to MNS in ipsilateral primary sensorimotor cortex (specifically Brodmann areas 3b and 1) are better explained by a GLM which takes account of natural variability in EEG responses than by a conventional GLM of consistent amplitude responses to MNS. Confirmation of this hypothesis would provide further evidence for a neural basis underlying the NBR, whilst identifying the features of the EEG response that correlate with natural single-trial NBR variability may provide further information about the functional significance of the NBR. In particular, since fluctuations in the power of alpha frequency oscillations are thought to reflect cycles of cortical inhibition (Klimesch et al., 2007;Mazaheri and Jensen, 2010), this work may improve the understanding of the relationship between inhibition and the NBR. We further hypothesize that natural trial-by-trial variability in the response to a consistent stimulus represents neuronal processes additional to those measured in the correlation between the average NBR magnitude and the intensity/ duration/frequency of the driving stimulus (Klingner et al., 2010;Shmuel et al., 2002).

Methods
Data were recorded on 18 right-handed subjects (age 27 ± 3 years, 8 female). This study was conducted with approval from the local ethics committee and informed consent was obtained from all subjects.

Paradigm
Stimulation was applied to the median nerve via two electrodes placed on the right wrist using square wave pulses of 0.5 ms duration (Digitimer DS7A, Letchworth Garden City, U.K.). The stimulation current amplitude was set (range 2.6-7 mA, mean 4.6 ± 1 mA) just above an individual's motor threshold so as to cause a small thumb distension. MNS was applied at 2 Hz in 40 blocks, each comprising a 10 s stimulation period followed by a passive rest period of 20.5-21 second duration. The onset of each MNS block was jittered by up to 500 ms to prevent stimulation always occurring at the same temporal latency relative to the MRI acquisition and to reduce the expectation of stimulus onset. Twenty individual MNS pulses were delivered in each stimulation block, giving a total of 800 pulses for each subject. The stimulation frequency of 2 Hz enabled distinct single-trial somatosensory evoked potentials (SEPs) to be recorded concurrently with ongoing oscillatory neuronal activity.

Data acquisition
EEG data were recorded using an MR-compatible EEG cap (EasyCap, Herrsching, Germany) with 63 electrodes following an extended international 10-20 system and an additional channel for recording the electrocardiogram (ECG). The reference electrode was positioned at FCz. A BrainAmp MR-plus EEG amplifier (Brain Products, Munich) with Vision Recorder (Version 1.10) was used for data acquisition. Data were recorded using a sampling rate of 5 kHz. The hardware band-pass filters were set to a 0.016-250 Hz range, with a roll-off of 30 dB/octave at high frequency. The impedance of all electrodes was kept below 20 kΩ. Gradient artefacts were minimised by mechanically isolating the EEG amplifiers from the scanner bed and minimising MR scanner room environment noise Mullinger et al., 2013a). In addition, the subject was positioned such that Fp1 and Fp2 were at the iso-centre in the foot/head direction so as to further reduce gradient artefact . Padding was placed around the subject's head to reduce motion-related artefacts. The EEG and MR scanner clocks were synchronised, and the TR made equal to a multiple of the EEG sampling period, to ensure consistent sampling of the waveforms (Mandelkow et al., 2006;Mullinger et al., 2008).
MR data were recorded in a Philips Achieva 3 T MR scanner (Philips Medical Systems, Best, Netherlands). Cardiac and respiratory cycles were simultaneously recorded using the scanner's physiological monitoring system (vector cardiogram (VCG) and respiratory belt). Firstly, a 5-minute simultaneous EEG-fMRI acquisition was used to localise S1/M1 and to check the EEG data quality (multi-slice EPI sequence: TR = 2 s, TE = 40 ms, 64 × 64 matrix, 3.25 × 3.25 mm 2 in-plane resolution, flip angle = 85°, SENSE factor = 2 and 20 slices of 5 mm thickness). Ten blocks of the MNS paradigm were delivered and IViewBOLD (Philips real-time GLM software) was used to calculate statistical maps of the BOLD responses. Based upon the location of the maximum BOLD response, the position of 10 contiguous axial slices covering S1/M1 for the subsequent main EEG-ASL-BOLD experiment was planned. A FAIR Double Acquisition Background Suppression (DABS) sequence (Wesolowski et al., 2009) was used for simultaneous acquisition of ASL and BOLD data (TR = 2.6 s, TE = 13 ms (ASL), 33 ms (BOLD), 2.65 × 2.65 × 5 mm 3 voxels, 212 mm FOV, SENSE factor = 2; post-label delay = 1400 ms, background suppression pulses at T BGS1 /T BGS2 = 340 ms/560 ms, scan duration~21 min). To facilitate co-registration and normalisation of functional data, a single volume EPI scan was acquired with the same geometry as the FAIR DABS images (TR = 8 s, TE = 40 ms, flip angle = 85°, SENSE factor = 2), along with both local, and whole-head anatomical image data-sets with 1 mm isotropic resolution. The locations of the EEG electrodes on the scalp surface were digitally recorded using a Polhemus isotrack 3D system (Polhemus, Vermont, USA) and co-registered with the subject's anatomical MRI scan.

Analysis EEG
Off-line EEG signal correction was based on averaging and then subtracting the gradient and pulse artefacts in Brain Vision Analyzer 2. A sliding-window gradient artefact correction was performed for each MR volume (5.2 s duration including tag and control image pairs) by subtracting an artefact template that was calculated by averaging over the 60 artefact repetitions temporally closest to the corrected volume. Cardiac cycles were defined from the VCG trace obtained from the MR scanner's physiological monitoring  and temporally aligned to the EEG data. A sliding-window pulse artefact correction was performed by subtracting an average template formed over 11 local cardiac cycles. Noisy channels and/or stimulation blocks were discounted by visual inspection, resulting in the exclusion of the data from three subjects from further analysis. Following artefact correction, data were down-sampled to 600 Hz and re-referenced to an average of all non-noisy channels. To allow separate quantification of the SEP and the ongoing oscillatory activity, the EEG signal was band-pass filtered into three frequency bands: 2-40 Hz, used to study SEPs; 8-13 Hz, used to study mu-activity (the mu rhythm is a typical example of the alpha rhythm, found in the sensorimotor cortex); and 14-30 Hz, used to study beta-activity. Data were exported to Matlab (MathWorks, USA) for further processing. As robust beta-band responses to MNS were not observed in all subjects, the subsequent analyses were focused upon SEP and mu data only.
A regularised, scalar beamformer was used to spatially localise brain responses to MNS (Van Veen et al., 1997), whilst also reducing sensitivity to residual gradient and pulse artefacts (Brookes et al., 2008(Brookes et al., , 2009. In order to localise the stimulus-related changes in both the SEP and mu power, pseudo-t-statistic (Ŧ) maps of SEP and mu responses to MNS were calculated over the whole head. The beamformer estimates the current source amplitude using a weighted sum of the voltage measurements at all EEG electrodes and a lead-field calculated from a triple-sphere head model for each location in the brain (0.5 cm isometric resolution) defined from each individual subject's anatomical MRI image. Subsequently the beamformer Ŧ-statistic is computed such that: where r denotes the position in the brain, P A represents the signal power during the 'active' stimulus-response window that contains the taskrelated response of interest, P P represents the signal power from a 'passive' window which defines a control or baseline period. N A and N P represent the noise power estimates inherent to the sensors during the active and passive periods respectively. For more details please refer to Brookes et al. (2008Brookes et al. ( , 2009).
SEP. the stimulus response window was defined as the 0.01-0.16 s timerange relative to each individual MNS pulse, with a passive window of 0.3-0.45 s (Brookes et al., 2009). Mu-power: the stimulus response window was defined as 0-9.5 s with a passive window of 20-29.5 s relative to the onset of each stimulus block. The mu-power passive window was chosen to ensure that effects of the post-stimulus rebound were not included in this window. Separately for SEP and mu-power Ŧ-maps, the location of the maximum change in activity between the stimulus response and passive windows in S1/M1 was chosen as the site of a virtual electrode (VE), from which a VE-time-course was extracted. At this stage, data from an additional two subjects were discounted from further analysis due to the presence of stimulus-linked motion artefacts. A further two subjects were discounted from the mu-power analysis due to the presence of large evoked-potential signals (identified by their time-and phase-locked nature) that corrupted the mu response. In total, 11 subjects (on average 37 ± 3 blocks) remained where all data quality was sufficient to perform the final analyses. The block-byblock changes in SEP and mu responses induced by MNS were quantified to facilitate comparison with fMRI signals. The temporal resolution of the fMRI signals was insufficient to allow for correlation with singletrial responses, therefore we focus on measuring the average response over each stimulus block.
SEP measurement. SEP VE-time-courses were epoched into single trials (− 50 ms pre-stimulus to 400 ms post-stimulus) and each trial was baseline corrected by subtracting the mean signal in the pre-stimulus period. The twenty baseline corrected SEP waveforms within each stimulation block were then averaged and the peak-to-peak P100-N140 amplitude of each block SEP waveform was measured using an automated linear regression method (Mayhew et al., 2006). Separately for each subject, the time-series of the 40 block SEP amplitudes was then mean-subtracted to form a regressor for subsequent GLM analysis. The amplitude of any rejected blocks was set to the mean value (zero).
Mu-power measurement. Mu VE-time-courses were Hilbert transformed and the average power in the stimulation (0-9.5 s) and passive (20-29.5 s) windows was calculated for each stimulation block. For each subject, a block mu-power regressor was then formed from the time-series of the 40 stimulation-period and 40 passive-period mupower amplitudes. Since our analysis focused on stimulus-induced changes in mu-power compared to the resting baseline, the mu regressor amplitude between 0 and 10 s was set to the mean amplitude of the active period and the mu regressor amplitude during the entire resting period (between 10 and 30 s) was set to the mean amplitude of the passive period (20-29.5 s). The amplitude of rejected blocks was set to the mean signal value.
fMRI pre-processing fMRI data were separated into BOLD and ASL data sets for subsequent analysis. The BOLD data were physiologically corrected using RETROICOR (Glover et al., 2000), whilst the use of background suppression meant that this correction was not required for the ASL data (Garica et al., 2005). All data were then motion corrected using FLIRT (Jenkinson et al., 2002) (FSL, http://www.fmrib.ox.ac.uk/fsl/) and linearly interpolated ("interp" function, Matlab, Mathworks USA) to an effective TR of 2.6 s. Tag-control ASL image pairs were then subtracted to create perfusion-weighted images. BOLD-weighted image pairs were averaged to produce mean BOLD-weighted data. Further pre-processing was carried out in SPM8 (http://www.fil.ion.ucl.ac.uk/spm/), where BOLD data were normalised to the standard MNI template and the same transform was applied to the ASL data. Data were smoothed with a Gaussian kernel (5 mm FWHM). A binary mask was created from the BOLD data to include only brain voxels that were present in all subjects. This was then used to mask each individual's BOLD and ASL data to ensure that GLM statistics were only calculated for voxels present in all subjects.
fMRI analysis GLM analysis of ASL data was carried out using a modified version of SPM5 (http://www.fil.ion.ucl.ac.uk/spm/). BOLD and ASL data were analysed separately. Data were modelled using a constant amplitude boxcar regressor of stimulation timings to encode the main effect of MNS, combined with a second regressor of either: 1) SEP block amplitude or 2) block mu power. All regressors were convolved with the canonical double-gamma SPM HRF. Both positive and negative T-contrasts were assessed for each regressor. Fixed-effects, second-level analyses were then performed to identify significant group-level activations separately for BOLD (p b 0.05, FWE corrected) and ASL (p b 0.001, uncorrected) data. Since the primary aim was to investigate S1/M1 responses, an MNI-space mask of S1/M1 cortex (Oxford-Harvard cortical atlas, FSL) was applied to all group level statistical maps. Subsequently, the grouplevel conjunction of the significant BOLD and CBF voxels for each contrast was used to mask individual subject's BOLD and CBF T-stat maps.

Region of interest (ROI) definition
Subject-specific, cubic regions-of-interest (ROI) (3 × 3 × 3 voxels) were centred on the peak BOLD T-stat voxel in S1/M1 for: (1) positive correlations with the boxcar (contralateral), (2) negative correlations with the boxcar (ipsilateral), (3) negative correlations with the SEP amplitude (ipsilateral) and (4) negative correlations with the continuous mu power (ipsilateral). BOLD and CBF, single-trial HR time-courses were then extracted from each of these BOLD ROIs, thereby allowing direct comparison between CBF and BOLD responses from the same spatial location. HR from individual trials were converted to percentage signal change relative to the final 6 s of each individual's mean HR across trials. For each ROI, BOLD and CBF HRs were then averaged across all subjects and the associated standard error across subjects calculated.

CMRO 2 estimation
The Davis model (Davis et al., 1998) is commonly used to estimate task-induced changes in CMRO 2 from concurrent measurements of the primary positive BOLD and CBF responses. Hypothesizing that the same physiological mechanism underlies both the PBR and NBR, we applied the Davis model to estimate the CMRO 2 change in both contralateral and ipsilateral S1/M1 cortex. For each subject, the peak signal change in the mean BOLD and CBF HR timecourses (extracted between 0 and 20 s) was found for both the positive and negative boxcar ROIs. The standard deviation of the signal over trials at the time point of the peak signal change was taken as the error in each of these measures. The change in CMRO 2 relative to baseline levels (CMRO 2 ) 0 was then calculated for each subject using Eq. (2): where the subscript 0 denotes signal values during the baseline period. α (the Grubb coefficient) was chosen to be 0.2, in line with recent MR literature (Chen and Pike, 2009), and β, which reflects deoxyhemoglobin concentration, was chosen to be 1.3 (Mark et al., 2011). M represents the maximum possible BOLD signal change: i.e. the change produced by an increase in CBF, which causes complete oxygen saturation in venous vessels. M is dependent on field strength and TE (Chiarelli et al., 2007) and is often calculated using a hypercapnic challenge (Gauthier et al., 2011). Since such a challenge was not performed in this study, a range of M-values appropriate for the sensorimotor cortex was taken from literature, normalised to field strength of 3 T and a TE of 33 ms, resulting in a value of M between 6% (Gauthier et al., 2011) and 14.9% (Kastrup et al., 2002). The percentage change in CMRO 2 (Δ%CMRO 2 ) was calculated for each of these M values as well as for the mean value of M (10.45%). For comparison with previous findings (Stefanovic et al., 2004) calculations with M = 9.5% (adjusted for field strength and TE) were also carried out. The calculated values of Δ%CMRO 2 were plotted against percentage change in CBF (Δ%CBF) for each subject and a linear fit was separately performed on the responses from the PBR and NBR regions to allow comparison of the values of the CMRO 2 /CBF coupling ratio (n) between the contralateral and ipsilateral regions.

EEG responses
Figs. 1A & B show the group average Ŧ-stat map of changes in both SEP (green) and oscillatory mu-power (purple) during the stimulus response window compared with the passive window. An increase in SEP power (denoted by positive Ŧ values) was only observed in contralateral S1 as expected from previous work (Hari et al., 1993;Kakigi, 1994) but a bilateral decrease in mu power (denoted by negative Ŧ values) was observed in response to MNS. The mu ERD during stimulation was found to be greater in the contralateral than the ipsilateral hemisphere, as previously reported (Yuan et al., 2010). Data from only two subjects showed local minima in the mu response in ipsilateral S1/M1, suggesting that the ipsilateral mu response exhibits much lower SNR than the contralateral S1/M1 mu response. Therefore to ensure consistent and robust response measurement we extracted responses from contralateral S1/M1. A high-degree of spatial overlap was observed between the SEP and mu responses in the contralateral hemisphere, although the group mean peak response location (shown by the crosshairs in Figs (Michel et al., 2004), it can be considered that these two signatures of neuronal activity arise from approximately the same cortical location. In each individual subject, the VE location of both SEP and mu responses were found in the contralateral hemisphere. The group average VE time-course of the SEP and mu responses to MNS from contralateral S1/M1 is shown in Figs. 1C & D, respectively. A clear SEP phase-locked to MNS (Fig. 1C) was observed with a positive peak at~100 ms and a negative peak at~180 ms. Mu power ERD was sustained during the entire stimulation period (0-10 s) (Fig. 1D), followed by recovery to pre-stimulus levels.

Haemodynamic correlates with the boxcar model
The T-stat maps of brain regions for which BOLD and CBF signals significantly correlated with the constant amplitude boxcar regressor formed using MNS timings are shown in Figs. 2Ai & Bi. We observed positive fMRI responses (red-yellow) to MNS in contralateral S1/M1 and negative fMRI responses (blue) to MNS in both ipsilateral S1/M1 and the supplementary motor area (SMA). A high degree of spatial correlation was observed between the regions of BOLD and CBF response. When assessing the spatial correlation of these responses it was important to ensure independence from differences in contrast-to-noise ratio (CNR) of BOLD and CBF data. Therefore the voxels with the top 5% of T-statistic values were taken separately for contralateral and ipsilateral S1/M1 for both BOLD and CBF responses (318 voxels for each region). The spatial overlap between the top 5% of BOLD and CBF voxels was found to be 64/52% for the PBR/NBR regions respectively. The group peak voxel locations (MNI co-ordinates, [x, y, z] mm) for the positive correlation with the boxcar model were [− 42 −22 46] for BOLD data and [−42 −16 42] for CBF data, whilst the group peak locations of negative correlation with the boxcar model were [34 −16 46] for BOLD data and [40 −18 46] for CBF data. The correlation with the constant amplitude boxcar model was more significant for the positive BOLD and CBF responses (peak T-stat values 18.8 and 11.1 respectively) than for the negative BOLD and CBF responses (peak T-stat values 12.4 and 6.9 respectively). Given that T-statistics reflect the goodness of the model fit to the data then these results suggest that either: the variability of the NBR amplitude over trials comprises a greater proportion of the overall signal in ipsilateral S1/M1 than the variability of the PBR in contralateral S1/M1, and/or the NBR amplitude is reduced relative to the PBR. It is also possible that a different canonical HRF may be required to optimally model the shape of the NBR compared to the PBR (Bagshaw et al., 2004). Fig. 3 shows the group averaged time-courses of BOLD and CBF responses to MNS, extracted from the positive and negative boxcar ROIs. Relative to baseline, we observe an increase in both the BOLD and CBF signals in the contralateral boxcar ROI and a decrease in the BOLD and CBF signals in the ipsilateral boxcar ROI. The BOLD and CBF response magnitudes are at least a factor of 2 greater in contralateral S1/M1 than in ipsilateral S1/M1. The ipsilateral BOLD and CBF responses exhibit longer peak latencies (10-12 s) than the contralateral responses (8-9 s).

Haemodynamic correlations with the SEP
Figs. 2Aii & Bii show the S1/M1 regions where the block-by-block modulations in SEP amplitude were negatively correlated with the BOLD and the CBF responses, respectively (overlaid in green). The group peak voxel location of BOLD-SEP correlation, ([36 − 20 46] (T = 8.8)), exhibited good spatial agreement with the region of CBF-SEP correlation, ([38 − 22 46] (T = 4.5)), in ipsilateral S1/M1. These SEP-fMRI correlations occurred in regions where negative fMRI responses to the boxcar model were also observed, as confirmed by the group mean time-courses extracted from the SEP ROIs (Fig. 3C). Therefore in this ipsilateral S1/M1 region we observe a negative correlation between SEP and the amplitude of the negative fMRI response. Alternatively, this can be conceptualised as a positive correlation between the SEP and the absolute magnitude of the ipsilateral S1/M1 fMRI response, meaning that stimulus blocks with highest SEP amplitudes are associated with the most negative (highest magnitude) NBR, see Supplementary Fig. S1 for a schematic illustration of this relationship. These negative responses displayed smaller peak magnitudes, but similar temporal profiles compared with the responses observed in the ipsilateral boxcar ROIs (Fig. 3B). No positive correlations between SEP amplitude and either the BOLD or CBF responses were observed in the S1/M1 network.

Haemodynamic correlations with mu power
Figs. 2Aiii & Biii show that mu power negatively correlated with block-by-block modulations in the BOLD and the CBF responses in ipsilateral S1/M1 and SMA regions of the S1/M1 (overlaid in purple). These mu-fMRI correlations also occurred in regions where negative BOLD and CBF responses to the boxcar model were observed, as demonstrated by the group average time-courses extracted from the mu ROIs (Fig. 3D). A schematic illustration designed to aid interpretation of the relationship between the NBR and the EEG responses is shown in Supplementary  Fig. S1. The magnitude of these negative responses was smaller than the responses observed in the ipsilateral boxcar ROI, but similar in magnitude to the response observed in the SEP ROI. The region of significant mu-BOLD correlation had a larger spatial extent than the region of significant SEP-BOLD correlation, and also had a slightly different peak location ([42 −12 46] (T = 7.5)). On further investigation, no significant correlation between the amplitude of the mu activity during stimulation and the SEP response was found in any subject (lowest p value = 0.11). A difference between the regions identified when using each of the two EEG regressors is therefore unsurprising. No positive correlations between mu power and either the CBF or BOLD responses to MNS were observed in the sensorimotor network.

CMRO 2 in PBR and NBR regions
The peak percentage signal change in the mean BOLD (Δ%BOLD) and CBF (Δ%CBF) response time-courses extracted from the positive (red) and negative (blue) boxcar ROIs are plotted for each subject in Fig. 4A. The black lines show iso-contours of constant CMRO 2 for M = 10.45%. These lines show that for most subjects an increase in CMRO 2 was associated with the PBR and a decrease in CMRO 2 associated with the NBR, this is also illustrated by Figs. 4B & C. We observe that the between-subject variability in BOLD and CBF data-points spans a range across several iso-contours of CMRO 2 which suggests that there was a change in the PBR region was not significantly changed if the outlying data point was removed (n = 0.5 ± 0.2). Fig. 4D displays the dependence of the coupling ratio (n) on the choice of the parameter M, and shows that n is always larger for the NBR region than for the PBR region, regardless of the M value used.

Discussion
Here, by using simultaneous EEG-BOLD-CBF recordings in humans we advance the understanding of the origin of the negative BOLD response and its relationship with concurrent measurements of cerebral oxygen metabolism and neuronal activity. We observed spatially coincident negative BOLD and CBF responses in S1/M1 regions of the cortical hemisphere ipsilateral to MNS and the amplitude of these ipsilateral S1/M1 responses was negatively correlated with the amplitude of both SEP and induced mu-power EEG responses ( Supplementary  Fig. S1). In ipsilateral S1/M1 the negative fMRI responses during MNS were accompanied by a decrease in CMRO 2 . There was also a significantly higher ΔCMRO 2 /ΔCBF coupling-ratio in ipsilateral S1/M1 than in the contralateral S1/M1 region, which exhibited a positive fMRI response to stimulation. These results suggest that an inhibitory neuronal mechanism underlies the NBR, but that the NBR is not simply the inverse of the PBR, since it exhibits a different neurovascular coupling.

Metabolic demand in PBR and NBR regions
We used the Davis model (Davis et al., 1998) to calculate the change in CMRO 2 evoked by MNS, thus allowing comparison of the metabolic demand in PBR and NBR regions. We found a significantly different value of the CMRO 2 /CBF coupling ratio, n, for regions with a NBR (0.91 ± 0.04) compared to those with a PBR (0.65 ± 0.03), assuming M = 10.45%. Whilst the values of n appear large compared with those previously reported (Shmuel et al., 2002;Stefanovic et al., 2004), this can be attributed to the parameter values employed here (α = 0.2 and β = 1.3), which are in line with most recent literature. If the "classical" value of α (α = 0.38, β = 1.3) is used with M = 9.5 (adjusted for our field-strength and TE) as previously derived for motor cortex (Stefanovic et al., 2004), then the resultant value of n = 0.48 ± 0.03 for the PBR region is in agreement with previous literature (Stefanovic et al., 2004). For these α and β values, the NBR region has a value of n = 0.79 ± 0.05, demonstrating that the difference in n between the PBR and NBR regions still remains for these "classical" values.
The difference in n that we observe between the PBR and NBR regions strongly suggests that a difference in the CMRO 2 /CBF coupling underlies the PBR and NBR. However, it is important to consider the sensitivity of the Davis model to the values of α, β and M when interpreting such differences in CMRO 2 /CBF coupling. CMRO 2 calculations have been shown to exhibit an increased sensitivity to the chosen value of α and β in NBR regions compared with PBR regions . Therefore the values employed in this study could result in an overestimation of CMRO 2 , and consequently n, in the NBR region. Theoretical values of α = 0.14 and β = 0.91 (M = 16.3%) have been derived which are designed to minimise the error in the calculated CMRO 2 for both PBR and NBR regions by abandoning the physical meaning of the parameters . Employing these theoretical values in this study resulted in n = 0.66 ± 0.03 for the PBR and n = 0.90 ± 0.04 for the NBR respectively. Since a significant difference between the coupling was still observed we suggest that the particular choice of α and β values does not cause the differences in n observed here. It is possible that there is a difference in CBF-CBV coupling between NBR and PBR regions, due to the different neuronal processes occurring in the two regions, meaning that different values of α should be used in the two regions. However, given the wide range of α-values tested here (0.14-0.38) and assuming that CBV and CBF remain coupled to some extent in both regions; it is unlikely that small differences in α between NBR and PBR regions could produce comparable n in the two regions.
The value of M employed also affects the calculation of n (Chiarelli et al., 2007), as shown in Fig. 4D. Calibrated BOLD measurements to estimate M for individual subjects were not performed in the current study. However, since both hemispheres are fed by similar vascular networks (Tatu et al., 1998) there is no reason to believe that the maximum BOLD response, and hence M, should vary significantly across hemispheres for S1/M1. Therefore, although it is not possible to definitively state the coupling ratio in the NBR and PBR regions, due to the dependence on the M value for S1/M1 cortex (Gauthier et al., 2011;Kastrup et al., 2002), Fig. 4D strongly suggests that the CMRO 2 /CBF coupling ratio n is significantly different for the NBR compared to the PBR region, regardless of the M value employed.
Previous studies found no significant difference in n between PBR and NBR regions (Shmuel et al., 2002;Stefanovic et al., 2004). Here, we find a value of n that is not significantly different to that obtained in previous studies (Stefanovic et al., 2004) if responses from the NBR and PBR regions are pooled (n = 0.64 ± 0.02). We suggest that the, previously unreported, significant differences which we observe between NBR and PBR regions may be due to either the differences in MR sequences used to collect the data (here a DABS pulse sequence was used at 3 T which allowed the simultaneous acquisition of BOLD and CBF data with high CNR for CBF and BOLD measures) or due to the MNS paradigm that we used.
Due to the stimulus duration employed, the fMRI responses measured in this study display peak signal changes within a short time window encompassed by 1 TR. Therefore we measure the peak amplitudes of the mean BOLD and CBF responses to MNS for each subject, rather than amplitudes averaged over a period of steady-state fMRI response which was used by previous studies (Kastrup et al., 2002;Stefanovic et al., 2004). However, given the number of blocks averaged over (37 ± 3) and the field strength (3 T) at which our data was acquired, there is no reason to believe that this analysis had a confounding effect on the experiment and the conclusions drawn.

Origin of altered oxygen metabolism-blood flow coupling in NBR regions
The difference in coupling between PBR and NBR regions is plausible given previous findings showing that n varies with attention (Moradi et al., 2012), caffeine administration  and stimulus duration (Lin et al., 2009). A variability in n is also consistent with the current hypothesis that changes in neuronal activity drive separate vascular and metabolic pathways that result in concurrent changes in CBF and CMRO 2 (Attwell and Iadecola, 2002;Buxton, 2010). This updates the more classical notion of a fixed mechanism whereby changes in CMRO 2 drive changes in CBF (Hoge et al., 1999;Shin, 2000). We hypothesise that the origin of the observed difference in n is related to the differences in the underlying neurophysiological basis of the PBR compared with NBR.
The CMRO 2 /CBF coupling is affected by the activity of excitatory/ inhibitory neurons and astrocytes with their associated energy demands and the release of signalling neurochemicals which regulate the dilation or constriction of the vasculature . However, the precise contribution of these factors and the interactions between them in vivo are currently unknown (Buxton, 2010;Lauritzen et al., 2012). The difference in the CMRO 2 /CBF coupling, that we observe, suggests that the NBR is not purely the inverse of the PBR (i.e. a simple reduction in excitatory neuronal firing), but rather caused by the recruitment of different classes of neurons/neuronal-astrocyte interactions. This could result in different metabolic (primarily driven by ATP cycles) and vascular (primarily driven by changes in intracellular Ca 2+ levels)  demands in the NBR region compared with the PBR region. This hypothesis is supported by evidence that unilateral sensorimotor stimulation induces NBR and also reduces the perception of threshold-level stimuli that are concurrently applied to the opposing hand (Klingner et al., 2010;Schafer et al., 2012). Taken together, the results presented here and in previous work suggest that NBR regions are associated with inhibitory neuronal processes that are not present in the PBR regions.

EEG-CBF response correlations
Close spatial agreement was observed in ipsilateral S1/M1 between the regions of negative EEG-BOLD correlations and negative EEG-CBF correlations. However, the clusters of CBF-SEP and CBF-mu correlations were both smaller in spatial extent and less statistically significant compared to the equivalent regions of EEG-BOLD correlation. The weaker correlation between EEG and CBF responses compared with EEG-BOLD responses may appear surprising, since it has been suggested that the arterial-capillary weighted CBF signal is more closely localised to neuronal activity than the venous BOLD signal (Duong et al., 2001), and so is potentially more sensitive to natural variations in neuronal activity. However, the lower CNR of CBF data compared to BOLD data likely explains the lack of sensitivity of CBF responses to the natural variations in the neuronal responses. The simultaneity of our recordings enabled us to rule out between-session differences as an explanation for the lack of tight coupling between EEG-CBF responses relative to EEG-BOLD responses, as previously suggested (Mayhew et al., 2010b). An alternative possible explanation is that since EEG is a macroscopic measurement of neuronal activity over a large network of neurons which reflects both the amplitude and synchrony of neuronal activity (Musall et al., 2012), it may be that trial-by-trial variability in EEG responses are more strongly correlated with the complete haemodynamic response, as represented by the BOLD signal.

Single-trial mu-fMRI correlations represent a link with cortical inhibition
Increases in the power of the mu oscillation (event related synchronisation (ERS)) and in the magnitude of the NBR (i.e. a more negative fMRI response) are both thought to reflect reduced cortical excitability or inhibition (Ferbert et al., 1992;Kastrup et al., 2008;Klimesch et al., 2007;Klingner et al., 2010;Liepert et al., 2001;Mazaheri and Jensen, 2010). In the current study we observe a negative correlation between the amplitude of the mu response and the amplitude of both BOLD and CBF responses in the NBR region, meaning that higher mu power is associated with increased magnitude of negative fMRI responses (Supplementary Fig. S1). Therefore, taken together with the changes in CMRO 2 in NBR regions discussed above, this trial-by-trial mu-fMRI correlation provides evidence that mu and negative fMRI responses possess a common link to the same inhibitory processes. Previous work in which it was demonstrated that stimuli delivered during high alpha power result in reduced visual PBR and enhanced NBR in the auditory cortex and default-mode network  is consistent with this hypothesis. Additionally, anti-correlations between the power of the 8-13 Hz posterior alpha oscillation and the resting-state BOLD signal have been widely reported in visual and parietal cortex (Goldman et al., 2002;Laufs et al., 2006) and are thought to represent cyclic fluctuations in the balance of cortical inhibition and excitation.
Previous studies in the sensorimotor system have reported taskdriven negative correlations between mu power and PBR, but no correlation with NBR (Ritter et al., 2009;Yuan et al., 2011). The negative mu-PBR correlation reported by Ritter et al. (2009) was primarily driven by the average reduction in mu power (ERD) and increase in BOLD signal in response to bilateral motor contraction, rather than natural response variability. Yuan et al. (2011) employed graded rates of unilateral finger tapping and observed a negative correlation between the contralateral PBR and both mu and beta responses (Yuan et al., 2011), which was again primarily driven by the scaling of the average BOLD response amplitude with the graded stimulus conditions. The apparent discrepancy between this previous work and the results of the current study are explained in the "Stimulus-driven modulations in the average response versus natural fluctuation in trial-by-trial amplitude" section below.
fMRI response modulations are widespread in the sensorimotor network Negative BOLD correlations with the boxcar model were observed both anterior and posterior to the contralateral PBR region (Fig. 2Ai, axial slice). A negative CBF response was observed anterior to the PBR (Fig. 2Bi, axial slice). The SMA also exhibited negative BOLD and CBF responses. Negative correlations between the amplitude of the NBR and the SEP were found in both the contralateral S1/M1 and SMA regions described above, whilst negative mu-NBR correlations were seen in the SMA. These results suggest that the trial-by-trial variability in NBR amplitude is similar across the whole sensorimotor network. We propose that the NBR regions of contralateral S1/M1 and SMA may represent a network of inhibitory activity associated with the primary ipsilateral S1/M1 NBR, sharing the same neurophysiological origin.
We believe that the trial-by-trial variability in the NBR which is modelled by the EEG regressors is related to neuronal effects during stimulation rather than slow haemodynamic responses due to the previous stimulus and the passive rest period of 20 s which was used. This is because the SEP regressor models only a measure of the neuronal response during the stimulation. Although the mu regressor does also model the variability during the passive resting period, if we model only the mu response during stimulation, we find very similar correlations to those presented in Fig. 2 for active and passive periods. Therefore, since EEG provides a direct measurement of neuronal activity with high temporal resolution, we believe it is unlikely that the correlations observed could arise spuriously due to the prior BOLD response.
The beamformer source localisation of the SEP response to MNS was strongly lateralised to contralateral S1/M1. The amplitude of the SEP extracted from this source was significantly negatively correlated with the negative BOLD response in the bilateral regions of sensorimotor cortex, meaning that the larger amplitude SEPs are associated with larger magnitude NBR. Upon further investigation, a negative SEP-PBR correlation was also observed in the PBR region of contralateral S1/M1 at a less significant statistical threshold (p b 0.001 uncorrected, data not shown). We therefore suggest that the boxcar model explains a greater proportion of the variance of the BOLD signal in contralateral S1/M1 than in ipsilateral S1/M1 because the amplitude of the PBR is larger and more consistent than that of the NBR (Fig. 3), which is supported by the Supplementary Tables S1 and S2. These results indicate that the SEP variability correlates with the BOLD response across the whole sensorimotor network, which is likely to be due to the underlying correlation between the trial-by-trial PBR and NBR amplitudes found in these data (Mayhew et al., 2013a).

Stimulus-driven modulations in the average response versus natural fluctuations in trial-by-trial amplitude
Previous work has shown that the average amplitude of: i) mu/beta power decreases relative to baseline (Muthukumaraswamy, 2010(Muthukumaraswamy, , 2011Stevenson et al., 2011;Yuan et al., 2011), ii) the SEP increases (Arthurs and Boniface, 2003), iii) the PBR increases relative to baseline Stevenson et al., 2011;Zhang et al., 2007) iv) the NBR decreases relative to baseline (Klingner et al., 2010;Shmuel et al., 2002) when the afferent stimulus input to the brain is driven with stimuli of increasing intensity, duration or frequency. These results may suggest that a positive SEP-PBR, negative mu-PBR (shown by (Yuan et al., 2011)) or positive mu-NBR correlation would be observed from studying average responses, compared to the negative trial-by-trial EEG-NBR correlations observed in this study. We hypothesise that the difference in the laterality of the BOLD (ipsilateral NBR, contralateral PBR) and EEG (bilateral mu ERD) responses and the negative EEG-fMRI correlations we observe may arise from different neuronal processes underlying extrinsic, stimulus-driven modulation of average response amplitude (Klingner et al., 2010) compared to intrinsic, natural trial-by-trial amplitude variability which contains fluctuations in subjects' internal processes such as attention and arousal (Macdonald et al., 2011;Prado et al., 2011;Sadaghiani et al., 2009Sadaghiani et al., , 2010. This may be supported by recent work showing that the PBR and the NBR amplitude are negatively correlated when measuring average response amplitudes driven by increasing stimulus intensity, but positively correlated when measuring natural amplitude variability in the trial-by-trial response to a consistent stimulus (Mayhew et al., 2013a).
The mu rhythm is not thought to be the primary driver of the metabolic demand that causes the BOLD signal changes during stimulation. Previous research suggests that neuronal activity across multiple frequency bands contributes to BOLD signal generation, with gamma frequency activity displaying the strongest correlation (Logothetis et al., 2001;Magri et al., 2012;Viswanathan and Freeman, 2007). Previous investigations of the superposition of EEG activity in different frequency bands have shown that the phase of the 8-13 Hz alpha oscillation modulates the power of high frequency (N 30 Hz) gamma activity (Osipova et al., 2008;Spaak et al., 2012). Moreover, correlations between alpha power and the BOLD response have been shown to occur during stimulation (Magri et al., 2012). Therefore we suggest that the variability in mu/SEP response amplitudes provides an index of changes in metabolic demand and neuronal activity, which are superimposed upon the primary source of the cortical response and therefore explain a substantial proportion of the variance in the fMRI response amplitude. This hypothesis of a superposition of activity, may also explain the lack of significant correlation of mu/SEP responses with the PBR. The PBR shows more significant correlation with the boxcar model than the NBR. Therefore the trial-by-trial variance of the PBR amplitude comprises a smaller proportion of the overall signal compared to the NBR, hence the PBR displays weaker correlation with the EEG responses than the NBR. This assertion is supported by the Supplementary Tables S1 and S2, which show the average percentage of variance explained by each regressor in each of the regions of interest.

NBR mechanisms
By using unilateral MNS we minimised the potential confound of blood-steal, since the ipsilateral NBR and contralateral PBR regions have different vascular territories. Our data suggest that a reduction in CBF with a lesser reduction in CMRO 2 is the most likely physiological origin of the NBR (Pasley et al., 2007;Shmuel et al., 2002;Stefanovic et al., 2004). However, our observation of differences in ΔCMRO 2 /ΔCBF coupling between contralateral and ipsilateral S1/M1, along with previous reports of delays in the onset and peak NBR latencies (Bagshaw et al., 2004;Klingner et al., 2011b), which are also seen here (Fig. 3), suggests that it may be too simplistic to explain the NBR by inverting the PBR neurovascular coupling mechanism. Furthermore, decreases in neuronal activity can arise through multiple mechanisms involving different contributions of decreased excitatory input and increased activity of inhibitory neuronal populations (Attwell et al., 2010;Boorman et al., 2010;Buzsáki et al., 2007;Cauli et al., 2004;Lauritzen et al., 2012;Logothetis, 2008;Schubert et al., 2008). Whilst inhibitory neuronal activity is potentially less energy demanding than excitatory activity (Buzsáki et al., 2007), recent work has demonstrated that it can in fact result in a PBR (Enager et al., 2009;Pelled et al., 2009). Further work is therefore required to fully understand the metabolic and vascular changes that accompany inhibitory activity. We believe that multimodal experimental approaches in humans that combine haemodynamic and neuronal measurements alongside mathematical modelling are essential to elucidate the origin of the NBR. In addition, valuable data could be obtained in animal models, by measuring fMRI responses simultaneously with separate recordings of the activity of local excitatory and inhibitory neuronal populations, in both PBR and NBR regions.
These additional studies would help to further explain the results presented here.
In conclusion, we report a negative correlation between NBR and both oscillatory mu power and SEP amplitude, providing a non-invasive method to relate EEG measures of neuronal activity and the NBR in humans. These correlations in conjunction with simultaneously measured reductions in CBF and CMRO 2 suggest that the NBR is, at least in part, neuronal in origin. However, differences in ΔCMRO 2 /ΔCBF coupling and time-course latency compared with the PBR data strongly suggest that the NBR does not simply originate from the inverse of the PBR mechanism, namely a reduction in excitatory activity. Therefore further detailed investigations are required to elucidate the precise mechanisms that underlie the NBR.