BOLD Consistently Matches Electrophysiology in Human Sensorimotor Cortex at Increasing Movement Rates: A Combined 7T fMRI and ECoG Study on Neurovascular Coupling

Abstract

Blood oxygenation level-dependent (BOLD) functional magnetic resonance imaging (fMRI) is widely used to measure human brain function and relies on the assumption that hemodynamic changes mirror the underlying neuronal activity. However, an often reported saturation of the BOLD response at high movement rates has led to the notion of a mismatch in neurovascular coupling. We combined BOLD fMRI at 7T and intracranial electrocorticography (ECoG) to assess the relationship between BOLD and neuronal population activity in human sensorimotor cortex using a motor task with increasing movement rates. Though linear models failed to predict BOLD responses from the task, the measured BOLD and ECoG responses from the same tissue were in good agreement. Electrocorticography explained almost 80% of the mismatch between measured- and model-predicted BOLD responses, indicating that in human sensorimotor cortex, a large portion of the BOLD nonlinearity with respect to behavior (movement rate) is well predicted by electrophysiology. The results further suggest that other reported examples of BOLD mismatch may be related to neuronal processes, rather than to neurovascular uncoupling.

Keywords

blood oxygenation level-dependent (BOLD) contrast electrocorticography (ECoG)fMRI neurovascular coupling sensorimotor cortex 7 tesla

INTRODUCTION

Neurovascular coupling is a key feature of the human brain. Models incorporating multiple electrical, cellular, chemical, and vascular parameters describe this relationship between the neuro-electrophysiological and the hemodynamic response to a stimulus in growing detail.^1–3 With these, abnormalities in coupling associated with brain disorders such as stroke and Alzheimer^2,4 may be understood better in terms of the underlying mechanisms. Accuracy of the models is obviously important, and therefore empirical data fed into the models need to be reliable. Functional magnetic resonance imaging (fMRI) using the blood oxygenation level-dependent (BOLD) signal is widely used to study human brain function as it can non-invasively detect and localize changes in brain state. Interpretation of fMRI is based on neurovascular coupling, as the BOLD signal represents the changing hemodynamics associated with neuronal activation. Deviation from a one-to-one relationship between neuronal activity and BOLD measures is a matter of concern, as it complicates interpretation of the models for brain pathology and of BOLD fMRI experiments.

There are several instances in which there is no simple one-to-one relationship between the BOLD signal, neuronal activity, and/or behavioral or stimulus properties.^5–9 An important example of this apparent mismatch is that the BOLD response amplitude saturates during rapid stimulus presentations or repeated movements at high movement rate.^10–12 The mismatch between motor execution and the BOLD response has fueled research on the underlying mechanisms elaborating on complex interactions between behavior, neuronal activity, and vascular signal sources such as blood flow, volume, and oxygen consumption.^1,13 For instance, a plausible source of saturation of the BOLD response is the ceiling effect where nonlinear BOLD signal behavior is caused by large changes in cerebral blood flow (CBF). A large increase in CBF is not accompanied by a proportionally large increase in BOLD signal change due to the washing out of all the deoxyhemoglobin in the venous vasculature. ^3,14 This mismatch is a matter of concern, as motor execution is a fundamental paradigm for many BOLD fMRI studies and is of clinical relevance in the prognosis of neurologic disorders such as Parkinson’s disease.^15,16 To date, the matter remains open to debate.

A principal challenge of previous work on relating BOLD responses to sensorimotor cortex activity is that, for the lack of direct measures of neuronal activity, behavioral measures are generally taken as a proxy. This assumes that there is a one-to-one relationship between sensory input and sensory cortex activity, as well as between movement rate and motor cortex activity. Recent findings however, showed that local neuronal activity measured with intracranial electrocorticography (ECoG) was suppressed with increases in hand movement rate, while the speed and amplitude of each hand movement remained the same.¹⁷ Thus, the assumption of a one-to-one relationship between movement rate and neuronal activity cannot be maintained. Given that at higher movement rates the BOLD response saturates, the finding that neuronal activity also declines raises the possibility that one may explain the other.

Here, we combined behavioral, BOLD fMRI, and direct measures of neuronal activity obtained with intracranial ECoG in the same human subjects, and examined whether BOLD responses accurately reflect neuronal activity over a range of different movement rates (0.3 to 2 Hz). Using ultrahigh-field fMRI (7 tesla), we obtained hemodynamic BOLD responses with high spatial and temporal fidelity and capillary specificity, weighted towards the hemodynamic demands of active neuronal tissue.¹⁸ Importantly, the ECoG electrode grids employed here allowed for a spatial resolution that is on the same order as our fMRI measurements. Studies have shown that power in higher frequencies (>40 Hz) correlates well with neuronal firing rates^19–21 and also with BOLD signal change.^22,23 The ECoG high-frequency broadband power thus provided a direct measure of electrical activity in confined neuronal ensembles, which we compared with BOLD fMRI measurements. We examined whether addition of the measured neuro-electrophysiological response in a model for the BOLD response affects the mismatch between model-predicted and – measured BOLD response. We computed expected BOLD time-courses for the movement (behavioral) and ECoG measures, by convolving both measures with an empirically obtained hemodynamic impulse response function (HRF). Two predicted hemodynamic BOLD time-courses were thus obtained, movement- and ECoG-based, which were compared directly with the measured BOLD time-course. The measured BOLD amplitude saturated at higher movement rates as expected. The movement-based BOLD amplitude revealed a substantial overestimation with respect to the measured BOLD amplitude at faster movement rates (≥ 1 Hz). The ECoG-based BOLD behaved quite similarly to the measured BOLD because of a considerable suppression in high frequency broadband power for these faster movement rates. The ECoG-based BOLD explained up to 80% of the discrepancy between the measured and the movement-based BOLD amplitude in somatosensory as well as motor cortex. Thus, even in the situation of considerable BOLD saturation and neuronal suppression, we observe a linear relationship between BOLD and ECoG HFB power over the whole range of executed movement rates in sensorimotor cortex.

MATERIALS AND METHODS

Subjects and Procedure

Three right-handed subjects (one woman, two men; mean age, 27 years; range, 19 to 43) participated in the study, approved by the Institutional Review Board of the Utrecht University Medical Center, after giving written informed consent in accordance with the Declaration of Helsinki 2008. The subjects had normal hand function and were scheduled for the implantation of ECoG arrays for the clinical purpose of epilepsy monitoring where the pathologic region did not extend to sensorimotor cortex. On average, 125 electrodes were implanted, 32 of which corresponded to a high-resolution grid of 8 × 4 electrodes located on the left sensorimotor cortex (Figure 1A). The high-resolution grid (pitch 3 mm) contained electrodes with 1 mm diameter, yielding a spherical measurement surface of 1 mm² per electrode. For each subject, the electrodes were localized using a method presented previously;²⁴ the electrodes were identified on a high resolution CT scan (Philips Tomoscan SR7000, Best, The Netherlands) and projected to the cortical surface rendering generated from a preoperative MRI scan on a 3T MR system (Philips 3T Achieva, Best, The Netherlands). All other functional and structural MRI data were obtained using a 7T MR System (Philips 7T Achieva, Cleveland, OH, USA).

Figure 1

Human sensorimotor cortex activation upon a hand movement task. Intracranial electrocorticography (ECoG) and 7T blood oxygenation level-dependent (BOLD) functional magnetic resonance imaging (fMRI) were combined to investigate the neuronal correlates of a BOLD saturation phenomenon, which is observed when increasing the hand movement rate (see panel C). (A) Blood oxygenation level-dependent activation of the sensorimotor cortex for the 0.3 Hz movement rate, together with the location of the high density ECoG grid on the 3D rendered brain of subject 1. Activation maps for the 0.3 Hz movement rate overlaid on the gradient-echo echo-planar imaging data is shown in Supplementary Figure S1A for all subjects together with the relative ECoG electrode positions (Supplementary Figure S1D). The primary motor (M1) and sensory (S1) cortices are also indicated. (B) From top to bottom for the 0.3-Hz movement rate; the obtained ECoG high-frequency broadband power (HFB power) responses for a single electrode, the hand movements as detected by the digital dataglove (from the ECoG session as shown in the top panel), and the corresponding BOLD response in primary motor cortex M1. Timing of the ECoG, dataglove, and BOLD data is with respect to the stimulus onset (t = 0 s). (C) The BOLD nonlinearity phenomenon in sensorimotor cortex for increasing movement rates. As previously reported, the BOLD amplitude saturates for higher movement rates (≥ 1 Hz.)²⁷ where it does not reach the amplitude as expected from a linear prediction model.

Functional Paradigm

The functional study was conducted in two sessions where the subjects performed the same tasks preoperatively during the 7T MRI session and postimplant during the ECoG session. In both sessions, the hand movements were recorded using a DataGlove 5 Ultra MRI (5DT, Irvine, CA USA, sampling rate 20 milliseconds), data of which were corrected for baseline drifts and converted to input functions for BOLD signal prediction analysis (see Data analysis paragraph below). The tasks consisted of closing the right hand on a visual cue from the rest position (open hand with palm facing up) at four different movement intervals: every 3, 2, 1, and 0.5 seconds, which translate to a movement rate of approximately 0.3, 0.5, 1, and 2 Hz (Figure 2A). The hand movements were executed during a 6-second period followed by a rest period of 15 seconds, and this is defined as one trial. In total, 15 trials were performed for each movement rate, where the order of movement rates was randomized across subjects.

Figure 2

Dataglove data illustrating the functional paradigm and the task performance for the blood oxygenation level-dependent (BOLD) functional magnetic resonance imaging (fMRI) and electrocorticography (ECoG) sessions. (A) Subjects were asked to close their right hand on a visual cue from the rest position (open hand with palm facing up) at four different movement rates: 0.3, 0.5, 1, and 2 Hz. The movements were performed during a 6-second period followed by a rest period of 15 seconds, which is defined as one trial. The hand movements were recorded using a dataglove during both the ECoG and fMRI session and gives a measure of the subjects’ performance for the different movement rates. Data shown are the hand movements as recorded by the dataglove for a representative trial for one subject during the BOLD fMRI session. (B) Movement rate as measured by the dataglove for both the BOLD (purple) and ECoG (red) session. Results shown are the mean ± standard deviation across all subjects. (C) Dataglove amplitude as measured for all movement rates for the BOLD (light green) and ECoG (dark green) session. The dataglove amplitude was similar across movement rates within a session. An amplitude difference was observed between the BOLD and ECoG sessions (consistent for all movement rates), potentially caused by different calibration of the dataglove on the subject’s hand. Results shown are the mean ± standard deviation across all subjects.

High-Field Imaging

Imaging was performed with a 16-channel head coil (Nova Medical, Wilmington, MA, USA) before surgical implantation of the ECoG arrays. All functional data were obtained using a multi-slice single-shot gradient-echo EPI technique with an isotropic voxel size of 1.5 mm, TR/TE = 880/27 milliseconds, flip angle = 65°, SENSE factor = 2.5, FOV = 155 × 155 mm², and 13 slices covering the left primary motor (M1), and primary sensory (S1) areas. To reduce image distortions, third order image-based shimming was performed (in-house developed IDL software, v6.3 RSI, Boulder, CO, USA) on the FOV of the functional scans after brain extraction. A high-resolution whole-brain T₂^*-weighted scan was acquired to visualize the large draining veins located on the pial surface, with the following parameters: 3D multishot gradient-echo echo-planar imaging, TR/TE = 90.8/27 milliseconds, flip angle = 20°, 2 averages, SENSE factor = 2, an isotropic voxel size of 0.5 mm, and FOV = 188 × 188 × 30 mm³. An EPI read out was used to have similar geometric distortions as in the functional images. This scan was used to manually create a draining vein mask, which was used to exclude the large draining veins located on the pial surface from the subsequent functional analysis (see Supplementary Figure S1C). The BOLD response of the large draining veins is known to be substantially delayed with respect to the BOLD response of the gray matter vasculature (the capillaries and intracortical veins and venules).^25,26 Including these, large draining veins could result in a less accurate, ‘blurred’, estimation of the ‘true’ gray matter BOLD response. A temporally biased gray matter BOLD response will affect the correspondence of the predicted BOLD signal in tracking the underlying neuronal activity. Each movement rate data set was registered to the movement rate data set that was acquired first, denoted as the ‘common space’, using Flirt (FSL, FMRIB Software Library, Oxford, UK). Activation maps were then obtained using FEAT: high-pass filtering at 1/21 Hz, no spatial smoothing, z-value threshold = 2 (cluster P threshold = 0.05) and convolution for the general linear model (GLM) using a standard double gamma variate hemodynamic response function (see Supplementary Figure S1A for the activation map for the 0.3 Hz movement rate). The area of intersection between the activation maps obtained for each movement rate produced a mask for each subject that contained voxels that were active in all movement rates. In addition, the areas M1 and S1 were manually delineated in the common space using anatomic boundaries, i.e., the central and pre and postcentral sulci. The intersection of the anatomic areas and the activation masks yielded the final M1 and S1 regions of interest (ROI) for further analysis (see Supplementary Figure S1B). Next, measured BOLD time-courses spanning all 15 trials were obtained for each movement rate after averaging across all voxels in the M1 and S1 ROIs followed by a ten-fold temporal interpolation. Finally, the BOLD time-courses were converted to units of percentage signal change by normalizing by a 3 s baseline period before the first trial.

Electrophysiology—ECoG Acquisition

ECoG data were acquired with a 128 channel recording system (Micromed, Treviso, Italy) with 512 Hz sampling rate and 0.15 to 134.4Hz band-pass filter and analyzed in the same way as in Hermes et al.¹⁷ Data were re-referenced to the common average of all electrodes recorded from the same amplifier. For subjects one and two, this computation included electrodes from the regular grids. To extract power changes in high frequencies, ECoG data were filtered for the 65 to 95 Hz bandwidth to obtain broadband gamma frequencies, using a third order Butterworth filter in two directions to minimize phase distortion (using the ‘filtfilt’ function in Matlab, The Mathworks Inc., Natick, MA, USA). After filtering, the log power of the analytic amplitude (by Hilbert transform) was calculated and the signal was smoothed with a 250-millisecond Gaussian window (standard deviation of 42 millisecond). The log power (HFB power) from the movement repetition tasks was normalized (z-score) with respect to the mean and standard deviation of all the 4-second period at the end of the trials. To map ECoG responses on M1 and S1, we selected the corresponding electrodes based on anatomic boundaries (the central sulcus and pre and postcentral sulci); electrodes located over a sulcus were not included. The z-scored log-power was then averaged across electrodes on M1 and S1. For all movement rates, we identified the individual ECoG responses within each trial and obtained the ratio between the first response and the rest of the responses for each trail. This ratio shows the neuronal suppression, if any, with movement rate. An average ratio per movement rate was obtained by averaging across trials.

Data Analysis

We modeled BOLD responses based on the measured movements (movement-based) and ECoG responses (ECoG-based) to determine whether addition of the measured neuro-electrophysiological response in a model can improve the prediction of the measured BOLD response. To investigate the BOLD signal linearity in relation to movement rate, we obtained the HRF for each subject from the slowest movement rate (0.3 Hz). As the electrophysiology data revealed very consistent neuronal responses for each movement at 0.3 Hz, we assume that the obtained HRF contains minimal detectable neuronal nonlinearities. We also assume minimal hemodynamic nonlinearities for the 0.3 Hz movement rate as previous studies observed a saturation of the BOLD response only at faster movement rates.^10,27 The 0.3 Hz HRF was estimated by deconvolution of the 0.3 Hz BOLD time-course with the corresponding dataglove input function using a conjugate gradient deconvolution method. The 0.3-Hz HRF was estimated separately for the M1 and S1 ROIs. The 0.3 Hz BOLD signal is thus modeled as:

S_{0.3 H z} = d_{0.3 H z} \otimes HRF + n

Where ⊗ denotes the convolution operator, S_0.3 _HZ the 0.3-Hz BOLD signal time-course, d_0.3 _HZ the 0.3-Hz dataglove trace, and n the noise. The HRF is then estimated by finding a least square solution in the Fourier domain using the conjugate gradient method:

HRF =_{HRF (f)}^{min} | | S_{0.3 H z} (f) - d_{0.3Hz} (f) \cdot HRF(f {)||}_{2}

where f denotes the Fourier domain.

Dataglove input functions were constructed from the postprocessed dataglove data using the timing and duration of the detected movements. The dataglove amplitude was not incorporated because a dataglove is potentially sensitive to calibration errors and different positioning on the subject’s hand, though it was consistent within the BOLD and ECoG sessions across movement rates (Figure 2). Convolving the HRF with the input functions from the dataglove data, we computed predicted BOLD time-courses for the other movement rates, referred to as ‘movement-based’ BOLD time-courses:

{\hat{S}}_{MB} = d_{MB} \otimes HRF

d_{MB} = \sum_{i}^{M} \sum_{j}^{N} A_{i, j} R_{i, j}

where A_i,j = 1, Ŝ_MB is the predicted ‘movement-based’ BOLD signal time-course for a particular movement rate (0.5, 1, or 2 Hz), d_MB the input function consisting of the amplitude A_i,j and the movement response shape Ryobtained from the dataglove data for corresponding movement rate. M is the number of trials, which was 15 for each movement rate, and N the number of movements within each trial.

Secondly, we modeled the drop in HFB power as a step function, expressed as the ratio between the peak amplitude of the first movement and the mean peak amplitude of the rest of the movements, averaged over all trials. The ECoG response ratios per movement rate obtained for each subject were used as weighting functions of the dataglove input functions. Convolving the HRF with these weighted input functions, we computed subject-specific ‘ECoG-based’ BOLD time-courses for each movement rate:

{\hat{S}}_{EB} = d_{EB} \otimes HRF

d_{EB} = \sum_{i}^{M} \sum_{j}^{N} A_{i, j} R_{i, j}

\begin{array}{l} w h e r e A_{i . j} = {\begin{cases} 1, j = 1 \\ r a t i o_{ECoG,} j \neq 1 \end{cases} \end{array}

Ŝ_EBis the predicted ‘ECoG-based’ BOLD signal time-course for a particular movement rate (0.5, 1, or 2 Hz), d_EB the input function where the A_i,j includes now the response ratios, ratio_ECOG, obtained from the ECoG session. The ECoG- and movement-based BOLD time-courses were compared with the measured BOLD time-courses for all movement rates where the goodness-of-fit was quantified in terms of the fit standard error (root mean square error, RMSE) and match in BOLD amplitude.

RESULTS

Figure 1 shows the experimental setup for the ECoG and BOLD measurements for one subject. Figure 1C illustrates the problem of investigation; the BOLD amplitude saturates for higher movement rates (> 1 Hz), and it does not reach the amplitude as expected from a linear prediction model. The dataglove data in Figure 2 shows that the subjects’ performance was similar between the movement rates within both the BOLD fMRI and ECoG sessions for all subjects. Figure 2A shows the dataglove traces during the BOLD fMRI session for one subject. Figures 2B and 2C shows the movement rate and movement amplitude as measured by the dataglove for both the BOLD fMRI and ECoG sessions for all subjects. For the BOLD session, the single ‘close-open’ hand movements were executed within 430 ± 24 milliseconds, 430 ± 90 milliseconds, 390 ± 90 milliseconds, and 300 ± 62 milliseconds for the 0.3, 0.5, 1, and 2 Hz movement rates respectively (mean ± standard deviation across subjects). For the ECoG session, the single ‘close-open’ hand movements were executed within 520 ± 220 milliseconds, 470 ± 160 milliseconds, 450 ± 145 milliseconds, and 280 ± 30 milliseconds for the 0.3, 0.5, 1, and 2 Hz movement rates, respectively (mean ± standard deviation across subjects).

BOLD Data Shows Saturation of the BOLD Response Amplitude for Movement Rates ≥ 1Hz

Average BOLD time-courses were obtained from voxels that were active for all movement rates for ROI in the primary motor (M1) and primary sensorimotor (S1) cortex for each subject (see Supplementary Figure S1B). For all subjects, the measured BOLD response amplitude saturated at movement rates ≥ 1 Hz, in contrast with both the predicted BOLD response based on the executed movements (Figure 3 for a representative subject) and the movements measured with the dataglove. The results of the measured BOLD response for all subjects for M1 and S1 are shown in Figure 6.

Figure 3

Blood oxygenation level-dependent (BOLD) amplitude shows a nonlinear relationship with movement rate. Measured BOLD responses (blue) saturate for movement rate ≥ 1 Hz with respect to the BOLD response predicted based on the executed movements (movement-based BOLD in red). Results shown for primary motor cortex (M1), averaged across trials, for a representative subject.

ECoG Broadband Gamma Power Shows Suppression of Amplitude for Movement Rates ≥ 1Hz

Measures of local neuronal activity (HFB power) were obtained for ECoG electrodes on the same ROIs in M1 and S1, as for the BOLD fMRI analysis for each subject. Figure 4A shows clear observable HFB power responses for each hand movement for a representative subject. For movement rates ≥ 1Hz and for all subjects we found a significant drop in HFB power for movements after the trial’s first movement as reported earlier.¹⁷ This indicates a suppression of neuronal activity with repeated movements at these rates, while the extent of measured movements remained the same (Figure 2C). Figure 4C shows the ECoG response ratios per movement rate obtained for all subjects for M1 and S1. For M1, the percentage drop in HFB power after the first movement was approximately 43 ± 14% for the movement rates ≥ 1 Hz. For S1, we found a percentage drop of approximately 35 ± 9% for the movement rates ≥ 1 Hz. For the movement rates < 1 Hz, we found a small reduction in HFB power, 12±9% and 9±10% of the HFB power of the first movement for M1 and S1 respectively (Figure 4C).

Figure 4

Electrophysiological data shows considerable suppression after the first movement for fast movement rates. (A) Electrocorticography (ECoG) data shown for electrodes on the primary motor cortex (M1) as measured by the log power in broadband gamma (HFB power) for each movement rate, 0.3, 0.5, 1, and 2 Hz, and averaged over trials for a representative subject (mean +/- standard). The time-courses show the average over the electrodes on M1. Z-scores are obtained by normalizing the data with the mean and standard deviation of the baseline period after each trial. Notice the drop in HFB power after the first movement for movement rates ≥ 1 Hz. Timing of the ECoG data is with respect to the stimulus onset (t = 0 seconds). (B) Location of the high-density ECoG grid shown in black, two electrodes from a regular ECoG are shown in gray, which were used to obtain HFB power responses from S1 for this subject. Electrocorticographic-grid locations for the other subjects can be found in Supplementary Figure S2. The dashed white line indicates the central sulcus. (C) The ratio of the first movement HFB power and all subsequent movements for all movement rates. Top and bottom panels show the results averaged over all electrodes located on primary sensorimotor cortex (M1 and S1) for all three subjects. Electrode locations were identified on a preoperative CT scan and were then projected to the cortical surface rendering generated from a preoperative magnetic resonance imaging (MRI) scan.²⁴

ECoG Responses Predict BOLD Responses in the Same Patch of Sensorimotor Cortex

Figure 5C shows how the dataglove data were used to obtain the movement-based and ECoG-based BOLD time-courses. When compared with the measured BOLD time-courses, we observed a large overestimation of the signal amplitude for movement rates ≥ 1 Hz for the movement-based BOLD (see red time-courses in the 1 and 2 Hz panels in Figure 5A for a representative subject). Supplementary Figure S2 shows the results for the other subjects. Figure 6 shows the mismatch of the movement-based BOLD with the measured BOLD amplitude and the increased fit standard error (RMSE) relative to the movement rate of 0.3 Hz, for all subjects. For the fast movement rates (≥ 1 Hz), we found that the movement-based BOLD amplitude overestimated the measured BOLD amplitude up to a factor of 1.9 and 1.8 for M1 and S1, respectively. For the 0.5-Hz movement rate, we found a factor of 1.2 overestimation of the measured BOLD amplitude using the movement-based BOLD model for both M1 and S1.

Figure 5

Prediction of blood oxygenation level-dependent (BOLD) time-courses using either a movement-based model or electrocorticography (ECoG)-weighted input functions. (A) BOLD primary motor cortex (M1) time-courses for all movement rates; measured (blue), as predicted from the movements (red, BOLD movement-based), as predicted by the subject-specific ECoG weighted input function (black, BOLD ECoG-based). Data is from the same subject as in Figures 1 and 4. Horizontal bars indicate the subject’s hand movements for each movement rate. Movement-based BOLD data was obtained by convolving the dataglove input functions with the hemodynamic impulse response function (HRF) estimated from the 0.3 Hz measured BOLD time-course. ECoG-based BOLD data was obtained by imposing the nonlinear neuronal suppression, approximated as a step function (computed ratios in Figure 4C), as a weighting function on the dataglove input followed by a convolution with the 0.3 Hz HRF. Time-courses for S1 for this subject are shown in Supplementary Figure S2. (B) Hemodynamic impulse response function as estimated from the 0.3 Hz BOLD time-course. (C) Illustration of how the behavioral responses as measured by the dataglove are transformed to neuronal responses and eventually the BOLD responses for the movement-based and ECoG-based approaches.

Figure 6

Blood oxygenation level-dependent (BOLD) amplitudes relative to the 0.3 Hz measured BOLD amplitude for the measured, movement- and electrocorticography (ECoG)-based BOLD time-courses in blue, red, and black respectively. M1 and S1 results on the left and right panel, respectively. ECoG-based BOLD shows a similar saturation behavior to the measured BOLD with respect to movement rate. Results averaged over all subjects for all movement rates. (A) Match of the predicted BOLD amplitude with the measured BOLD amplitude, mean +/- standard deviation across all subjects; measured (blue), movement-based BOLD (red), and BOLD amplitude as predicted by the subject-specific ECoG weighted input function (black). The movement-based BOLD overestimates the measured BOLD amplitude up to a factor of 1.9 for the faster movement rates (≥ 1 Hz). The ECoG-based BOLD reduced the movement-based BOLD overestimation by approximately 82 ± 3% for M1 for all movement rates ≥ 0.5 Hz. For primary sensory cortex (S1), the electrocorticography (ECoG)-based BOLD reduced the movement-based BOLD overestimation by 50%, 78%, and 70% for movement rates 0.5, 1 and 2 Hz, respectively. (B) Goodness-of-fit of the predicted BOLD time-courses to the measured BOLD expressed as the fit standard error root mean square error (RMSE). RMSE of the ECoG-based BOLD are significantly reduced compared with the movement-based RMSE for movement rates ≥ 0.5 Hz in M1 (P<0.03) and for movement rates ≥ 1 Hz in S1 (P<0.002, one-sided paired Student’s t-test).

Figure 5A (black time-courses) shows for a representative subject the ECoG-based BOLD time-courses, utilizing the electrophysiological response ratios from the subject’s ECoG session (Figure 4C) as weighting functions. Supplementary Figure S1 shows the results for the other subjects. We found that the ECoG-based BOLD closely matched the measured BOLD time-courses in both the M1 and S1 ROIs for each subject (see Figure 6). Including the ECoG weighted input functions into the model reduced the overestimation of the modeled BOLD amplitude by 82 ± 3% for M1 as compared with the movement-based BOLD for all movement rates ≥ 0.5 Hz. For S1, the ECoG-based BOLD reduced the movement-based BOLD overestimation by 50%, 78%, and 70% for movement rates 0.5, 1, and 2 Hz, respectively. The accuracy of the ECoG-based approach in predicting the measured BOLD is also reflected by the lower RMSE values; the RMSE values of the ECoG-based approach were significantly reduced with respect to the RMSE values of the movement-based approach for movement rates ≥0.5 Hz in M1 (P<0.03) and for movement rates ≥ 1 Hz in S1 (P<0.002, one-sided paired Student’s t-test).

DISCUSSION

We obtained and combined intracranial ECoG, 7T BOLD fMRI, and behavioral measures from the same human subjects for a simple motor task with several rates of movement. Goal was to evaluate the relationship between BOLD and neuronal population activity in sensorimotor cortex, in the context of new evidence suggesting that overt movements are not tightly linked to cortical activity.¹⁷ ECoG and movement measures were used to construct BOLD prediction models, which were compared with the measured BOLD time-courses. The saturation of the BOLD response with respect to movement rate as reported in literature^{10,11,27–29} was confirmed in the present study. Importantly, the prediction model based on the HFB power in the ECoG proved to yield a BOLD time-course that was strikingly similar to the measured BOLD time-course. This indicates that the BOLD response is tightly coupled to neuronal activity over the range of executed movement rates. This finding directly impacts on the validity of inferring neuronal activity patterns in sensorimotor studies from BOLD-fMRI data. Moreover, this finding suggests that observed BOLD nonlinearities in other cortical regions or experimental conditions may also relate to the underlying neuronal processes rather than neurovascular uncoupling.

In this study, BOLD signals were acquired with high spatial detail from sensorimotor cortex using a 7T scanner. ECoG offers well-localized electrophysiological measurements of neuronal activity, and electrophysiological data were obtained without the use of anesthetics. The latter, often applied in animal studies, poses a potential contamination of the neurovascular coupling mechanism³⁰ and thus adds uncertainty in the obtained measures. Simultaneous measurements of 7T BOLD fMRI and ECoG were not possible for safety reasons (risk of electrode heating due to radio frequency energy deposition), but the use of a digital dataglove made it possible to bridge the two modalities. The recently revealed nonlinearity between hand movements and the underlying neuronal activity¹⁷ led us to reconsider the assumption that at increasing movement rates, neuronal responses remain invariant and that therefore saturation of the BOLD response in sensorimotor cortex is due to neurovascular uncoupling. The ECoG HFB power shows considerable decline in neuronal activity after the initial movement for the movement rates ≥1 Hz for both M1 and S1.¹⁷ As broadband HFB power changes have been linked to neuronal population firing rate,^19–21 the observed decline in the ECoG HFB power at higher movement rates in this study could suggest that, after movement initiation, fewer motor cortex neurons have to fire or firing rates may possibly decline when multiple similar movements are made at a fast rate.¹⁷ Including this subject-specific drop in neuronal activity as a weighting function in a BOLD fMRI prediction model, we obtained a close match with the measured BOLD data from the same cortical regions in the same subjects (Figures 6A and 6B).

The fact that a linear relationship between BOLD and ECoG HFB power was found in both M1 and S1 suggests that the principle may generalize to other brain regions and functions also. Previous studies that have combined electrophysiology and BOLD or blood flow measurements have focused mainly on sensory cortices such as the visual or somatosensory cortex using passive neuronal stimulation, either invasively in rodents or primates, or non-invasively using scalp-based electrophysiology measurements electroencephalography (EEG) in humans.^5,31–37 While nonlinearities in hemodynamic and electrophysiological measurements are often reported, two studies on human visual cortex using modulation of the stimulation frequency suggested that the hemodynamic response nonlinearity can potentially be explained by the nonlinear electrophysiological activity as captured by the EEG event-related potentials.^31,32 Electroencephalography measurements, however, do not allow for accurate localization of the source of electrophysiological responses because of the inverse problem.³⁸ Moreover, event-related potentials such as those measured by EEG do not map well onto BOLD activity for changing stimulus duration as was recently shown with ECoG,³⁹ whereas changes in HFB power do, as reported in the same and other data sets.^17,40 Electrocorticography measures directly from active neurons confined to cortical patches of similar spatial extent as BOLD fMRI voxels, with a close correlation between HFB power and BOLD signal change data sets.^17,40 Hence, the combination of ECoG and BOLD fMRI in humans has more straightforward and direct implications for relating BOLD to underlying neuronal activity.

While we found a tight relationship between BOLD and HFB power for M1, this seems to be slightly less tight in S1. For the movement rate of 2 Hz, the ECoG-based BOLD could only explain 70% of the mismatch between the movement-based and measured BOLD amplitude in S1, whereas this was 85% in M1 (Figure 6). This remaining saturation in the measured BOLD amplitude in S1 could be explained by neuronal processes that are not captured by the HFB power but by lower frequency bands.^23,41Although we found a significant decrease in the HFB power for subsequent movements at high movement rates, we observed no movement rate effect for the amplitude in the lower frequency, beta-band power, response (12 to 28 Hz).¹⁷ The beta-band activity, however, remained suppressed between movements at faster movement rates in both motor and sensory cortex, in contrast with the slower movement rates where the beta-band activity recovered to baseline values.¹⁷ Also, the underlying neuronal activity in sensory and motor cortex differs during motor action. The HFB response is a single clear transient increase in power in M1, whereas in S1, it consists of two sequential increases associated with closing and opening of the hand.¹⁷ Also, sensory cortex activity may be affected by presynaptic inhibition of afferent input, which is not likely to occur in the motor cortex. This inhibition may change with movement rate. Alternatively, BOLD nonlinearities that are caused by differences in the local vasculature and the creation of the BOLD signal^26,42 could still contribute.

The observed BOLD nonlinearities could potentially be caused by the ceiling effect where the nonlinear saturation is due to large changes in cerebral blood flow (CBF) that wash out the deoxyhemoglobin resulting in a plateauing of the BOLD signal amplitude.^13,14 Earlier PET studies, however, have observed similar saturation of CBF changes in sensorimotor cortex for fast movement rates (up to 2 Hz), as we observed for the BOLD signal.^27,28 Additionally, a PET study by Dettmers et al⁴³ compared BOLD changes with CBF changes in motor cortex and observed a linear relationship between BOLD and CBF in a finger tapping task (1 Hz) for a range of tapping forces.⁴³ These studies indicate that it is very likely that our observed BOLD changes are linearly related to the underlying blood flow changes, and thus contain a very minimal ‘ceiling effect.’ For future work, it would be interesting to obtain CBF data for example provided by arterial spin labeling MRI techniques to validate previous findings. Nevertheless, the present study indicates that for the sensorimotor system, the impact of these vascular processes on the BOLD signals is much less than proposed in literature. The present findings also encourage further investigation to assess the linearity between BOLD and neuronal activity for other brain regions and stimulus configurations, such as continuous cortical activation or ranges of stimulus intensities and durations. Similar investigations would be interesting for animal studies that are typically performed to investigate neurovascular coupling.

In summary, in the present study, the relationship between BOLD and neuronal population activity in the human sensorimotor cortex was explored over a range of movement rates, given new evidence of a mismatch between behavioral measures and electrophysiology. Here we investigated relationships between behavioral, hemodynamic, and direct cortical electrophysiological measures, in a human sensorimotor experiment. Our data show considerable saturation of the BOLD response with respect to hand movement rate and neuronal response suppression after the first movement, in the presence of rate-invariant movement measures. When including the subject-specific suppression in ECoG HFB power in our BOLD fMRI prediction model, we obtain a high goodness-of-fit of our measured BOLD data for all movement rates. This provides evidence that a large portion of the BOLD saturation is well explained by the drop in HFB power, and thus the neuronal processing associated with the HFB power. Thus, while there can be a discrepancy between executed movements and both the BOLD and ECoG responses, the latter are closely matched within the same patch of sensorimotor cortex. Further research is warranted to translate these finding to other cortices and brain functions.

DISCLOSURE/CONFLICT OF INTEREST

The authors declare no conflict of interest.

Footnotes

ACKNOWLEDGEMENTS

We thank Frans Leijten, Cyrille Ferrier, Peter van Rijen, and Peter Gosselaar for enabling data acquisition in the patients, and Mariska van Steensel, Erik Aarnoutse, and Martin Bleichner for ECoG data collection. We also thank David Leopold for his helpful suggestions on the manuscript.

Notes

Supplementary Information accompanies the paper on the Journal of Cerebral Blood Flow & Metabolism website ()

References

Lauritzen

Mathiesen

Schaefer

Thomsen

. Neuronal inhibition and excitation, and the dichotomic control of brain hemodynamic and oxygen responses. Neuroimage 2012; 62: 1040–1050.

Attwell

Buchan

Charpak

Lauritzen

Macvicar

Newman

. Glial and neuronal control of brain blood flow. Nature 2010; 468: 232–243.

Cardoso

Sirotin

Lima

Glushenkova

Das

. The neuroimaging signal is a linear sum of neurally distinct stimulus- and task-related components. Nat Neurosci 2012; 15: 1298–1306.

Girouard

Iadecola

. Neurovascular coupling in the normal brain and in hypertension, stroke, and Alzheimer disease. J Appi Physiol 2006; 100: 328–335.

Devor

Dunn

Andermann

Ulbert

Boas

Dale

. Coupling of total hemoglobin concentration, oxygenation, and neural activity in rat somatosensory cortex. Neuron 2003; 39: 353–359.

Soltysik

Peck

White

Crosson

Briggs

. Comparison of hemodynamic response nonlinearity across primary cortical areas. Neuroimage 2004; 22: 1117–1127.

Vazquez

Noll

. Nonlinear aspects of the BOLD response in functional MRI. Neuroimage 1998; 7: 108–118.

Ances

Zarahn

Greenberg

Detre

. Coupling of neural activation to blood flow in the somatosensory cortex of rats is time-intensity separable, but not linear. J Cereb Blood Flow Metab 2000; 20: 921–930.

Sirotin

Das

. Anticipatory haemodynamic signals in sensory cortex not predicted by local neuronal activity. Nature 2009; 457: 475–479.

10.

Jancke

Specht

Mirzazade

Loose

Himmelbach

Lutz

. A parametric analysis of the ‘rate effect’ in the sensorimotor cortex: a functional magnetic resonance imaging analysis in human subjects. Neurosci Lett 1998; 252: 37–40.

11.

Miezin

Maccotta

Ollinger

Petersen

Buckner

. Characterizing the hemodynamic response: effects of presentation rate, sampling procedure, and the possibility of ordering brain activity based on relative timing. Neuroimage 2000; 11(6 Pt 1): 735–759.

12.

Sadato

Ibanez

Campbell

Deiber

Le Bihan

Hallett

. Frequency-dependent changes of regional cerebral blood flow during finger movements: functional MRI compared to PET. J Cereb Blood Flow Metab 1997; 17: 670–679.

13.

Buxton

Uludag

Dubowitz

Liu

. Modeling the hemodynamic response to brain activation. Neuroimage 2004; 23(Suppl 1): S220–S233.

14.

Miller

Luh

Liu

Martinez

Obata

Wong

. Nonlinear temporal dynamics of the cerebral blood flow response. Hum Brain Mapp 2001; 13: 1–12.

15.

Ebersbach

Hattig

Schelosky

Wissel

Poewe

. Perseverative motor behaviour in Parkinson’s disease. Neuropsychologia 1994; 32: 799–804.

16.

Redgrave

Rodriguez

Smith

Rodriguez-Oroz

Lehericy

Bergman

. Goal-directed and habitual control in the basal ganglia: implications for Parkinson’s disease. Nat Rev Neurosci 2010; 11: 760–772.

17.

Hermes

Siero

Aarnoutse

Leijten

Petridou

Ramsey

. Dissociation between neuronal activity in sensorimotor cortex and hand movement revealed as a function of movement rate. J Neurosci 2012; 32: 9736–9744.

18.

Yacoub

Shmuel

Pfeuffer

Van De Moortele

Adriany

Andersen

. Imaging brain function in humans at 7 Tesla. Magn Reson Med 2001; 45: 588–594.

19.

Ray

Maunsell

. Different origins of gamma rhythm and high-gamma activity in macaque visual cortex. PLoS Biol 2011; 9:e1000610.

20.

Miller

Sorensen

Ojemann

den Nijs

. Power-law scaling in the brain surface electric potential. PLoS Comput Biol 2009; 5:e1000609.

21.

Manning

Jacobs

Fried

Kahana

. Broadband shifts in local field potential power spectra are correlated with single-neuron spiking in humans. J Neurosci 2009; 29: 13613–13620.

22.

Lachaux

Fonlupt

Kahane

Minotti

Hoffmann

Bertrand

. Relationship between task-related gamma oscillations and BOLD signal: new insights from combined fMRI and intracranial EEG. Hum Brain Mapp 2007; 28: 1368–1375.

23.

Hermes

Miller

Vansteensel

Aarnoutse

Leijten

Ramsey

. Neurophysiologic correlates of fMRI in human motor cortex. Hum Brain Mapp 2012; 33: 1689–1699.

24.

Hermes

Miller

Noordmans

Vansteensel

Ramsey

. Automated electrocorticographic electrode localization on individually rendered brain surfaces. J Neurosci Methods 2010; 185: 293–298.

25.

Lee

Glover

Meyer

. Discrimination of large venous vessels in time-course spiral blood-oxygen-level-dependent magnetic-resonance functional neuroimaging. Magn Reson Med 1995; 33: 745–754.

26.

Siero

Petridou

Hoogduin

Luijten

Ramsey

. Cortical depth-dependent temporal dynamics of the BOLD response in the human brain. J Cereb Blood Flow Metab 2011; 31: 1999–2008.

27.

Blinkenberg

Bonde

Holm

Svarer

Andersen

Paulson

. Rate dependence of regional cerebral activation during performance of a repetitive motor task: a PET study. J Cereb Blood Flow Metab 1996; 16: 794–803.

28.

Sadato

Ibanez

Deiber

Campbell

Leonardo

Hallett

. Frequency-dependent changes of regional cerebral blood flow during finger movements. J Cereb Blood Flow Metab 1996; 16: 23–33.

29.

Riecker

Wildgruber

Mathiak

Grodd

Ackermann

. Parametric analysis of rate-dependent hemodynamic response functions of cortical and subcortical brain structures during auditorily cued finger tapping: a fMRI study. Neuroimage 18: 731–739.

30.

Martin

Martindale

Berwick

Mayhew

. Investigating neural-hemodynamic coupling and the hemodynamic response function in the awake rat. Neuroimage 2006; 32: 33–48.

31.

Zhang

Liu

Chen

. Noninvasive study of neurovascular coupling during graded neuronal suppression. J Cereb Blood Flow Metab 2008; 28: 280–290.

32.

Wan

Riera

Iwata

Takahashi

Wakabayashi

Kawashima

. The neural basis of the hemodynamic response nonlinearity in human primary visual cortex: implications for neurovascular coupling mechanism. Neuroimage 2006; 32: 616–625.

33.

Ureshi

Matsuura

Kanno

. Stimulus frequency dependence of the linear relationship between local cerebral blood flow and field potential evoked by activation of rat somatosensory cortex. Neurosci Res 2004; 48: 147–153.

34.

Rees

Friston

Koch

. A direct quantitative relationship between the functional properties of human and macaque V5. Nat Neurosci 2000; 3: 716–723.

35.

Logothetis

Pauls

Augath

Trinath

Oeltermann

. Neurophysiological investigation of the basis of the fMRI signal. Nature 2001; 412: 150–157.

36.

Janz

Heinrich

Kornmayer

Bach

Hennig

. Coupling of neural activity and BOLD fMRI response: new insights by combination of fMRI and VEP experiments in transition from single events to continuous stimulation. Magn Reson Med 2001; 46: 482–486.

37.

Sheth

Nemoto

Guiou

Walker

Pouratian

Toga

. Linear and nonlinear relationships between neuronal activity, oxygen metabolism, and hemodynamic responses. Neuron 2004; 42: 347–355.

38.

Grech

Cassar

Muscat

Camilleri

Fabri

Zervakis

. Review on solving the inverse problem in EEG source analysis. J Neuroeng Rehabil 2008; 5:25.

39.

Engell

Huettel

McCarthy

. The fMRI BOLD signal tracks electrophysiological spectral perturbations, not event-related potentials. Neuroimage 2012; 59: 2600–2606.

40.

Huettel

McKeown

Song

Hart

Spencer

. Linking hemodynamic and electrophysiological measures of brain activity: evidence from functional MRI and intracranial field potentials. Cereb Cortex 2004; 14: 165–173.

41.

Scheeringa

Fries

Petersson

Oostenveld

Grothe

Norris

. Neuronal dynamics underlying high- and low-frequency EEG oscillations contribute independently to the human BOLD signal. Neuron 2011; 69: 572–583.

42.

Duvernoy

Delon

Vannson

. Cortical blood vessels of the human brain. Brain Res Bull 1981; 7: 519–579.

43.

Dettmers

Connelly

Stephan

Turner

Friston

Frackowiak

. Quantitative comparison of functional magnetic resonance imaging with positron emission tomography using a force-related paradigm. Neuroimage 1996; 4: 201–209.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

1.67 MB

0.02 MB