Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Neuroimage. Author manuscript; available in PMC 2010 April 1.
Published in final edited form as:
PMCID: PMC2848457

Study of neurovascular coupling in humans via simultaneous magnetoencephalography and diffuse optical imaging acquisition


By combining diffuse optical imaging (DOI) and magnetoencephalography (MEG) we investigate neurovascular coupling non-invasively in human subjects using median-nerve stimulation. Previous fMRI studies have shown a habituation effect in the hemodynamic blood oxygen level-dependent (BOLD) response for stimulation periods longer than 2 s. With DOI and MEG we can test whether this effect in hemodynamic response can be accounted for by a habituation effect in the neural response. Our experimental results show that the habituation effect in the hemodynamic response is stronger than that in the earliest cortical neural response (N20). Using a linear convolution model to predict hemodynamic responses we found that including late neural components (≥30 ms) improves the prediction of the hemoglobin response. This finding suggests that in addition to the initial evoked-response deflections related to the talamic afferent input, later cortical activity is needed to predict the hemodynamic response.

Keywords: Magnetoencephalography (MEG), Diffuse optical imaging (DOI), Near infrared spectroscopy (NIRS), Neurovascular coupling, Somatosensory response, Primary sensory area (SI)


Several functional magnetic resonance imaging (fMRI) studies have shown a strong non-linear relationship of the vascular response to variations in stimulus duration (Vazquez and Noll, 1998; Birn et al., 2001; Birn and Bandettini, 2005). In particular, Birn et al. (2001), Bandettini et al. (2002), and Birn and Bandettini (2005) have shown that, while the hemodynamic responses to long-duration stimuli behave in an approximately linear fashion, the responses to short-duration stimuli (<2 s) are larger than those predicted from a linear model. The non-linearity of the blood oxygen level-dependent (BOLD) signal could arise from the habituation effect in the neuronal response, the habituation effect in the hemodynamics, or both. Since the above fMRI experiments did not involve simultaneous electrophysiological measurements, the origin of the non-linear vascular responses remained unclear.

To date, a large effort has been devoted to studying the relationship between the neural and the vascular evoked responses in animals measuring the two types of signals simultaneously with invasive techniques (Ngai et al., 1999; Jones et al., 2001; Norup Nielsen and Lauritzen, 2001; Devor et al., 2003; Sheth et al., 2003; Logothetis and Wandell, 2004; Ureshi et al., 2004; Martindale et al., 2005; Royl et al., 2006). Non-invasive studies of the neurovascular coupling are more challenging in humans because of the difficulty in simultaneous acquisition of BOLD fMRI and electroencephalography (EEG) data. While progress has been recently made (Negishi et al., 2004; Ritter and Villringer, 2006; Schmid et al., 2006; Vasios et al., 2006; Mantini et al., 2007; Riera et al., 2007), the signal-to-noise ratio (SNR) of the combined measurements is still not optimal. The use of near-infrared spectroscopy (NIRS)/diffuse optical imaging (DOI), instead of fMRI to measure the vascular response, in combination with EEG/magneto-encephalography (MEG) achieves an excellent SNR since there is no interference between the optical and electrical measurements. NIRS/DOI has recently been introduced to study neurovascular coupling in humans in combination with EEG (Obrig et al., 2002; Butti et al., 2006; Rovati et al., 2007; Koch et al., 2008; Herrmann et al., 2008) and DC-MEG (Mackert et al., 2008).

In this work, we developed a simultaneous MEG-DOI recording technique to examine the neurovascular relationship using somato-sensory stimulus trains with varying durations. With DOI, we observed a non-linear vascular response similar to that described in Boynton et al. (1996), Birn et al. (2001), and Birn and Bandettini (2005): while the hemodynamic responses to 2-s or longer stimulus trains behave in an approximately linear fashion, the response to a 1-s train is larger than that predicted from a linear model. To test whether this vascular non-linearity for a short stimulus train can be explained by the habituation effect in the neuronal response, we acquired MEG data simultaneously with DOI. These two modalities are compatible with each other, and they allow for simultaneous measurement of neural and vascular responses with a few technical precautions. Despite its lower spatial resolution, DOI provides a higher temporal resolution and good SNR for small stimuli such as the median-nerve stimulus compared to fMRI.

We chose median-nerve electrical stimulation because the corresponding somatosensory responses are thoroughly studied in intracranial electrophysiological experiments in animals. Therefore, our results can be directly compared with, for example, the many forepaw stimulation studies in rats. In addition, the electrical response to median-nerve stimulation has been extensively studied in humans with EEG and MEG(Brenner et al.,1978; Goff et al.,1978; Kaufman et al.,1981; Hari et al., 1984; Tiihonen et al., 1989; Allison, 1992; Mauguiere et al., 1997). It has been shown that a median-nerve stimulus activates a complex somatosensory cortical network (Hari and Forss, 1999), including the contralateral primary sensory cortex (cSI), the bilateral secondary sensory cortices (SII), and other associated areas such as the posterior parietal cortex (PPC). Unlike SII, which lies in the inferior parietal region, deep in the highly folded cortex, the SI hand area, located on the posterior wall of the central sulcus, is relatively superficial and is easily studied with both MEG and DOI. PPC, located on the wall of the post-central sulcus, medial and posterior to the SI hand area, is also relatively superficial to the cortex compared with SII. With DOI we recover both oxy- and deoxy-hemoglobin (HbO and HbR) evoked responses in cSI for different stimulus train durations. From MEG data, we estimate the sites of neural activity and extract the cSI current dipole time course. By co-registering the two measurements with respect to each subject's structural MRI, we verified that there is spatial co-localization between neural and vascular activations.

The MEG response to median-nerve stimulation typically contains four main peaks: N20, P35, N45, and P60. The first peak N20, occurs at 20 ms after stimulus onset, corresponds to intra-cellular currents towards the surface of the cortex, induced by the early excitatory postsynaptic potentials (EPSPs) (Wikström et al., 1996). The cSI response then changes its polarity and reaches a maximum at 35 ms (P35), corresponding to intra-cellular currents flowing inward from superficial layers into the cortex. This peak has been suggested to be due to the early inhibitory postsynaptic potentials (IPSPs). P35 is followed by a narrow negative peak N45, corresponding to secondary EPSPs, and a broader positive peak P60. N20 and P60 are known to be robust with respect to stimulus repetition rate. In contrast, P35 and N45 are highly sensitive to stimulus repetition rate. For instance, several studies have reported that N45 is enhanced with increase in the stimulus repetition rate (Wikström et al., 1996). In our case, using a stimulus repetition rate of 4 Hz, we observed a strong habituation of these two peaks during the stimulation trains.

The high SNR of optical and MEG measurements allows us to compare different MEG peaks, and evaluate which best predicts the hemodynamic responses as measured by DOI. By limiting the consideration to specific peaks of the neural responses, we aim to isolate the part of neuronal activity that causes the local excitatory hemodynamic response. We compare the linear predictions of the hemoglobin responses using stimulus onsets, the current dipole amplitude of a specific peak, the sum of the peak amplitudes, and the root mean square (RMS) value of the current dipole amplitude in a time window before the response returns to the baseline.

This paper is organized as follows: Methods section describes our simultaneous MEG-DOI acquisition technique and mathematical models in the data analysis. Results section presents our experimental results in the spatial and the temporal domains, followed by a discussion.


Experimental details

We studied 6 healthy, right-handed subjects (5 males, 1 female). The subjects were 24–61 years old (mean 37 ± 14). Informed consent was obtained from subjects in accordance with the Massachusetts General Hospital ethical committee prior to participation. Measurements were conducted in a magnetically shielded room, on subjects laying on a non-magnetic bed. During the measurements, subjects were alone in the room with the lights dimmed.

Trains of median-nerve stimuli (slightly above motor threshold, current <10 mA, 0.2-ms pulse duration, 4-Hz repetition rate) were delivered at the right wrist according to an event-related protocol, with a random inter stimulus interval (ISI) varying between 2 and 18 s (mean 12 s). The train durations were 1, 2, 3, and 4 s, corresponding to 4, 8, 12, and 16 stimuli, respectively. These four conditions were presented in a pseudo-random order (Dale, 1999) in each of the six runs (6 minutes per run), with a minimum of 2-minute rest period between consecutive runs depending on the subject's condition. Ideally, we would like to explore a larger range of stimulus train durations; however, we have observed that the MEG response to stimulus trains beyond 4 s is often contaminated by eye-movements and blinks.

The DOI data were acquired with a multi-channel continuous-wave optical imager (CW4 system, TechEn Inc.) as described in Franceschini et al. (2003). The CW4 instrument was kept in the control room and we used 10-m long fibers to deliver and collect light to and from the subjects in the magnetically shielded room. The source and detector fiber bundles were attached to three black rubber stripes and arranged in two rows of eight detectors separated by a row of nine sources, for a total of 32 source–detector combinations with 3-cm separation constituting the probe array (see Fig. 1a). In each of the nine sources, two laser diodes emitted light at 690 nm and 830 nm wavelength, respectively. The laser wavelengths were chosen to maximize the SNR of the measurements, while minimizing the cross-talk error between HbO and HbR (Yamashita et al., 2001; Strangman et al., 2003; Sato et al., 2004; Uludag et al., 2004). The sources were frequency-encoded by steps of approximately 200 Hz between 4.0 kHz and 7.4 kHz, and their signals were acquired simultaneously by the sixteen avalanche photodiode detectors. The output from each detector was digitized at 40 kHz. The demodulation of the source signals was performed off-line by using an infinite-impulse-response filter, and finally the signals were converted to 10 Hz time series for analysis.

Fig. 1
DOI probe. (a) Planar view of the DOI source/detector arrangement. (b) The DOI probe mounted on a subject's head. (c) Schematics of a custom-made plug (1 cm thickness) that houses one fiber optic bundle.

The optical probe covered an area of 6 × 16 cm2 and was positioned to extend from left frontal to the left parietal areas with the center approximately above Brodmann area c3 (Fig. 1b). The fiber tips terminated at 90° bend and were mounted to the probe by means of low-profile custom-made plugs (1 cm thickness) in order to fit the probe in the limited space between the head of the subject and the MEG dewar (Fig. 1c). All materials used in the optical probe were tested in the MEG environment, and special fiber jackets and epoxy, different from those typically used for optical alone or the optical-MRI measurements, were used.

The MEG data were acquired simultaneously with the optical data using a 306-channel Neuromag VectorView MEG device (Elekta-Neuromag Oy, Helsinki, Finland), in a magnetically shielded room (Imedco AG, Switzerland) equipped with adjustable lighting, an intercom, and a video camera. The VectorView sensor array consists of 204 planar gradiometers and 102 magnetometer detectors. Eye movements and blinks were monitored with electrooculogram (EOG). The MEG signals were recorded continuously with a 1-kHz sampling rate and minimally filtered (0.01 to 250 Hz). There was no special setting adjustment needed for MEG recording in a simultaneous MEG-DOI acquisition, but the SNR of MEG signals suffers slightly from the additional distance between the sensors and the head introduced by the optical probe. The positions of the optical probes and the four head-position indicator coils for MEG were digitized before the measurement with a 3D magnetic digitizer (Polhemus).

In a separate session, MRI anatomical images were collected with a Siemens Avanto 1.5 T scanner with the T1-weighted sagittal MPRAGE protocol and the fast low angle shot (FLASH) protocol. We segmented the head surface and the gray matter cortex from the T1-weighted images (Dale et al., 1999; Fischl et al., 1999). The head surface was then co-registered with the digitized points obtained in the joint MEG-DOI session. From the FLASH images, we extracted the inner skull surface to construct the MEG forward matrix using a single-compartment boundary-element model (Hämäläinen and Sarvas, 1989; Oostendorp and van Oosterom, 1989).


The optical data were temporally co-registered with the MEG signals by collecting the stimulation triggers with the two systems. The optical raw data were processed off-line using HomER (, an in-house software with a graphical user interface for visualization and analysis of NIRS data implemented in MATLAB (Mathworks Inc., Sherborn, MA). Data were band-pass filtered between 0.02 and 0.5 Hz to reduce the effect of slow drifts and heartbeat, and the temporal changes in the intensity were translated into temporal changes in absorption coefficients using the modified Beer-Lambert law. Concentrations of HbO and HbR were derived from the absorption changes at the two wavelengths using hemoglobin extinction coefficients reported in (Wray et al., 1988). To reduce motion artifacts and systemic oscillation, principal component analysis (PCA) was applied by removing 1–3 components, as described in Zhang et al. (2005). Finally, the hemodynamic responses to each stimulus train condition were calculated using an ordinary least squares linear deconvolution. The deconvolved hemoglobin responses were averaged over the measured runs to improve the SNR with respect to single-run responses.

Due to the sparse distribution of our source–detector pairs, we limited our analysis to a simple back-projection reconstruction method to identify the vascular activation area. From the 32 source–detector combinations, we identified the source–detector pair with the most statistically significant deoxy-hemoglobin activation (p-value for the selected source–detector pairs across six subjects <0.05), and used the data of this combination to investigate the neurovascular coupling relation. For the chosen source–detector pair, we computed the hemoglobin time courses for each subject and each condition. Compared to HbO, HbR is known to be less sensitive to systemic oscillations, such as respiration and arterial pulsation (Obrig et al., 2000; Franceschini et al., 2003). We performed identical analysis for HbR and HbO; however, due to the larger contamination to HbO, and hence a worse SNR, we only display the figures related to HbR in Results section.

For MEG data, we first performed baseline correction and spatial whitening based on a noise covariance matrix computed from signals in a 200-ms window before the stimulus onset. With the MEG forward model and source space constructed, we can estimate the neural response from the MEG measurements by solving the inverse problem (Baillet et al., 2001). In this work, we employed both dipole fitting (Mosher et al., 1992) and a standard distributed inverse solver, the minimum-norm estimate (MNE) (Hämäläinen and Ilmoniemi, 1984), constrained to the cortical gray matter (Dale and Sereno, 1993; Dale et al., 2000). The standard dipole fitting algorithms estimate the location, orientation, and amplitudes of a fixed number of current dipoles. Dipole fitting is known to be robust for source signal estimation; however, localization is challenging when several sources are active. On the other hand, the MNE solution is the smallest-energy source distribution in the source space that reproduces the MEG measurements. Our study of the neurovascular coupling in the spatial domain is based on the MNE results. The analysis in the temporal domain is based on the estimated time courses from a standard dipole fitting procedure (Nelder and Mead, 1965; xfit software).

In dipole fitting, we separately fitted the evoked responses of each of the four duration conditions, using a selection of MEG channels (about 45 channels) covering an area above the left SI cortex, at 20–38 ms after the onset of the first stimulus in the train. Subsequently, the time course of this source over the whole epoch was determined with a linear least squares fit. In our single-dipole model, the estimated SI dipole time course usually includes contribution from nearby regions such as PPC. It may also include later activations which often have a more distributed spatial pattern. We believe the neurovascular coupling analysis is still valid despite the fact that the estimates may contain both SI and PPC activations. That is because the hemodynamic optical signals often reflect a larger spatial pattern, which may very likely contain integrated effects including SI and PPC.

For each subject, after averaging runs for a particular duration condition, we identified N20, P35, N45, and P60 train response as indicated in Fig. 4 for condition one in subject 1. While several studies have reported that N45 strengthens with reduced stimulus repetition rate, for the 250-ms setting in our experiments, N45 is relatively weak and its amplitude is further affected by P35 due to temporal proximity. Hence, we omitted this peak from the remaining analysis.

Fig. 4
Subject S1's estimated neural responses to the 1-s stimulus train duration condition (four stimuli). The stimulus onsets and the detected peaks are marked. For the neurovascular coupling analysis, we considered each individual peak amplitude, the difference ...

Neurovascular coupling analysis

To confirm the non-linearity of the vascular responses when the stimulus duration was varied as shown in previous fMRI studies (Birn et al., 2001; Bandettini et al., 2002; Birn and Bandettini, 2005), we applied a linear model to the optical data to predict the responses to long stimulus train duration using the response to short stimulus train. For instance, to linearly predict the response to a 2-s stimulus train, we added the response to a 1-s train and its 1-s delayed version. To linearly predict the response to 3-s stimulus train, the 1-s response was replicated three times, shifted by 0, 1, and 2 s, and then summed. We performed similar linear predictions for the response to a 4-s train using the 1-s and the 2-s responses.

To evaluate the correlation between the hemoglobin and neural response, we followed two data analysis approaches. In the first approach, we evaluated the linearity of the responses by plotting the response amplitudes versus the number of stimuli and calculated the coefficient of determination (R2) to compare results. In particular, for the hemodynamic response and for each subject, we calculated the area under the HbR and HbO responses for each condition and normalized them with respect to the largest response area. This normalization removes artifacts due to subject variability in head size and in probe location with respect to center of activation. For each condition, we then averaged the responses of the six subjects together and evaluated linearity with zero crossing imposed, as well as the zero intersection without imposing zero crossing, for HbO and HbR, as shown in Fig. 5. For the MEG responses, for each subject and each condition, we determined the amplitudes of N20, P35, and P60 responses to each stimulus in the train (i.e., amplitude responses to 4 stimuli for condition one, 8 stimuli for condition two, and so on) by calculating the maximum (or minimum) in a narrow time window. Since the signals were baseline-corrected prior to dipole fitting, the amplitudes of the peaks were set relative to the zero baseline. To remove artifacts due to subject variability in source depth and head size, we normalized the peaks amplitudes with respect to the amplitude of the largest P35 in response to the first stimulus in the train among the four conditions. We then averaged the cumulative MEG responses to each stimulus across subjects, deriving the first 4 responses from condition one, the second 4 responses from condition two, responses to stimuli 9 to 12 from condition three, and finally responses to stimuli 13 to 16 from condition four, and we calculated R2 and zero intersection of the three peaks, N20, P35, and P60.

Fig. 5
Average normalized responses of HbR (area) and average normalized amplitudes of the three peaks in the neural responses across subjects under the four stimulus duration conditions.

In addition to the peak amplitudes, we also considered the sum of the peak amplitudes in absolute value, including P35 − N20, P60 − N20, and P35 + P60 − N20, and the root mean square (RMS) value of the dipole amplitude in the time window after N20 and before the response returning to baseline1 (see Fig. 4). The chosen time window covers the early activations such as SI and PPC (Hari and Forss, 1999). We also performed calculations including different time windows and did not see statistical differences on the results with respect to the chosen one.

In the second more rigorous approach, for each subject we estimated the linear convolution model between the measured hemoglobin time series to each condition and the measured neural response time series, including N20, P35, P60, sum of peak amplitudes in absolution values, and the RMS value, as described above.

For all experiments, the temporal length of the impulse response function (h) was set to 10 s, which corresponded to 100 time samples for the hemoglobin responses. We estimate a separate impulse response function for HbO (hHbO(t)) and HbR (hHbR(t)). This convolution model, also used in Logothetis and Wandell (2004) and Martindale et al. (2005), assumes that the neuronal activity and hemoglobin responses are related by


where N(t) is the time course of neuronal activity represented by a MEG component for each stimulus pulse in the pulse train, yHbO(t) and yHbR(t) are the time courses of the hemoglobin signals,2 and hHbO(t) and hHbR(t) are the hemoglobin response functions. In the above, [large star] indicates the convolution of two time series. nHbO(t) and nHbO(t) include background physiology and measurement noise, which we assume to be mean zero Gaussian white noise. The hHbO(t) and hHbO(t) are estimated simultaneously for all stimulus conditions using least squares linear deconvolution for each subject separately, and we denote the corresponding estimates as h HbO (t) and h HbO (t). Predicting the variation in HbO and HbR responses to the different stimulus conditions is then possible using a linear transformation. Results for each subject are shown in Fig. 6.

Fig. 6
Linear convolution analysis. For the four stimulus train duration conditions and for each of the six subjects, the measured HbR using DOI (black curve) and the predicted HbR based on the linear convolution model with different components in the neural ...

A good linear fit is reflected by a good match between the predicted and measured hemoglobin response, which is mathematically evaluated using the coefficients of determination (R2). For example, in the analysis of HbR, the corresponding R2 becomes


We then convert the R2 value into the Fisher's z-score:


A two-sample t-test based on the Fisher's z-scores between the linear predictions from two types of neural components were calculated to assess the statistical significance (p-values ≤ 0.05 between different components).


Co-localization of MEG Sources and DOI Signals

The median-nerve stimulus employed in this work is known to robustly activate cSI (Brodmann area 3) (Hari and Forss, 1999). To verify the co-localization of the neural and the vascular activations, we present the reconstructed activation maps from MEG data (yellow area) using the MNE software (Hämäläinen, 2005) and the contour levels of HbR signal changes overlaid on the corresponding MRI cortices (Fig. 2). The reconstructed HbR map for each subject was obtained in HomER using a standard back-projection method. The contour levels represent the 90%, 75%, and 50% thresholds of maximum signal decrease between 3 and 5 s post-stimulation onset with respect to −2-0 s before stimulation onset for the 4-s stimulus train duration condition. The optode positions recovered with the 3D digitizer are depicted as red and blue circles in Fig. 2. For each subject, the source–detector pair with the most statistically significant HbR changes are shown in blue. In one subject (S6), the digitization of the optical probe was lost due to repositioning of the probe after digitization, and only the MEG source location is shown. In the other five subjects with digitized optode positions, the vascular responses, reflected by the largest signal decreases in HbR, are spatially co-localized with the estimated neural activations from the MEG measurements. The results for HbO, not shown, are consistent with the ones for HbR.

Fig. 2
Activation maps for six subjects. Neural activation estimated from MEG measurements at 35 ms after the onset of the first stimulus in the stimulus train is depicted in yellow. The contour plot corresponds to 90%, 75%, and 50% thresholds of maximum signal ...

Non-linear DOI responses

The grand-averaged HbRs over all six subjects obtained from the selected source–detector pairs for 2-, 3-, and 4-s train durations are presented in Fig. 3 (solid curves). We used the HbR responses corresponding to 1- and 2-s trains separately to linearly predict the responses to longer stimulus train duration as described in the previous section. The linear predictions (dashed curves) are overlayed with the corresponding measured HbR. We can see that the HbR response to a 1-s train overestimates the responses to longer train durations using a linear model, reflecting an overshoot of vascular supply to short stimulus trains. In contrast, the HbR response to a 2-s train provides a good linear fit to the response to the 4-s train. Our DOI results confirm that the responses to short stimulus trains are larger than those estimated using a linear model (Bandettini et al., 2002; Birn et al., 2001; Birn and Bandettini, 2005).

Fig. 3
Linear predicted HbR from measured HbR to short stimulus duration. Solid curves show the grand-average HbR measured with DOI, and dashed curves are the corresponding predictions assuming a linear model. Top: linear prediction of responses to 2-, 3-, and ...

Linear MEG responses

Fig. 4 shows the evoked cSI response to the 1-s stimulus train in a subject (S1), with the major peaks highlighted. The evoked responses to stimulus 5 through stimulus 16 in the other duration conditions are similar to the response to the fourth stimulus in the 1-s condition, but slightly noisier.

For most of the subjects, the amplitudes of N20 and P60 are similar in all stimuli in a train. On the other hand, P35's amplitude drops significantly in response to the second to the last stimuli in a train. If the vascular response were linearly related to either N20 or P60, the vascular response would have been linearly related to the stimulus train duration for the four conditions that we examined. In contrast, if the vascular response was linearly related to P35, the vascular response would have a similar amplitude across the four conditions. Therefore, the linear relationship between HbR and any single peak is not adequate.

Fig. 5 shows the cumulative sum of the amplitudes of each of the three peaks averaged over the 6 subjects for each train duration. In contrast to the hemodynamic responses (also shown in Fig. 5 as the grand average of the areas under the HbR responses for each condition) which is not linear with respect to train duration (R2 = 0.96), both N20 and P60 are linear with respect to stimulus train duration (R2 ≥ 0.99). On the other hand, P35 exhibits habituation effect in response to the second and later stimuli in a train, reflected by the smallest R2 value (R2 = 0.93) and the largest zero intersection (0.26).

Hemodynamic predictions

Fig. 6 and the Supplementary materials illustrate the measured HbR, yHbR(t), and the linear prediction, h HbR (t)[large star]N(t), using peak amplitudes, sum of the peak amplitudes in absolute value, and the RMS value of the neural response for the six subjects and the four stimulus train durations. It reveals that none of the peaks can individually explain the HbR response with a linear convolution model. It also reveals that the sum of three peaks (red dashed curves) provides relatively better linear fit to HbR consistently across all six subjects.

We evaluated the quality of the predictions based on differences in the coefficients of determination (R2) between the predicted and the measured time courses of the hemodynamic responses. Table 1 summarizes the results of a two-sample t-test. In particular, Table 1 (top) shows a cross-reference of p-values for statistically significant differences in R2 values corresponding to different components in predicting HbR. Table 1 (bottom) shows the corresponding results for HbO. The average R2 across six subjects is presented in the right-most column. The sum of peak amplitudes and the RMS value of the neural response provide the best linear prediction compared to any single component, reflected by the largest R2 values. All of them are statistically better than the input stimuli in predicting HbR, and the sum of the three peaks also achieves statistically better prediction than N20.

Table 1
Cross-reference of p-values for significant differences in R2 values for each subject and each component used to predict HbR (top) and HbO (bottom).

Since HbO is more sensitive to systemic physiological effects than HbR, the SNR is lower than that for HbR. Although the results for HbO are similar to those for HbR, the R2 values of the HbO predictions are smaller than the corresponding ones for HbR prediction, and statistically significant better prediction is established for the sum of three peaks and the RMS value with respect to N20.


Previous studies using MEG and NIRS in combination (Mackert et al., 2004; Sander et al., 2007; Mackert et al., 2008) have employed direct current (DC) MEG measurement techniques to study DC or near-DC changes in the magnetic fields in a block design finger-movement paradigm with 30-s blocks of finger movements alternating with rest periods. Limited to one stimulus type (stimulus frequency and duration), these studies cannot explore the linear/nonlinear neurovascular coupling relationship. In contrast, here we investigate the relationship of the NIRS/DOI vascular signals with the detailed time course of MEG current sources. The use of an event-related rather than a block design allows us to employ shorter stimulation periods (1–4 s versus 30 s) and to evaluate the contribution of single neural component to the hemodynamic responses. Our approach enables us to directly compare the results with the vast number of neurovascular coupling findings in animals using invasive electrical and vascular recordings.

The choice of limiting the study to four conditions is dictated by the fact that in order to examine more duration conditions, it is necessary to extend the scanning time to ensure a sufficient SNR in the average signals for the neurovascular coupling analysis. The choice of 1 s as the shortest condition is due to the fact that hemoglobin response for shorter stimuli is too noisy; the choice of 4 s as the longest duration is due to the fact that MEG signal for longer stimuli become noisier because of eye blinks and muscle contractions. In any case, the 1 to 4 s duration chosen is within the range of non-linearity previously observed in the hemodynamic response with BOLD fMRI (Bandettini et al., 2002).

We observed a good spatial agreement between the activation areas determined by the two imaging modalities. The neural source location determined using MNE in five subjects is in good proximity with the maximum hemoglobin evoked changes measured with NIRS. This result reveals the benefit of integrating other vascular measurements such as fMRI, which has a good spatial resolution, to help solve the ill-posed MEG inverse problem. Our result also indicates that a combined MEG-DOI inverse formulation may provide additional information for the two ill-posed inverse problems.

We verified that the hemodynamic response is non-linear with respect to stimulus train duration (R2 = 0.96 and 0.22 for zero intersection). Our result shows that the hemodynamic responses to short stimulus trains (<2 s) measured with NIRS are larger than those predicted by a linear model, and it is in agreement with fMRI results in humans (Bandettini and Ungerleider, 2001; Birn et al., 2001). With the simultaneous MEG measurements, we tested whether the hemodynamic non-linearity can be justified by a habituation effect in the neural responses. For the tested stimulus train durations, we found that N20 and P60 show a linear relationship with train duration, achieving R2 value above 0.99 with 0.26 zero intersection. On the other hand, due the strong habituation effect in P35, the R2 value for P35 is 0.93; however, this habituation effect is too strong to explain the habituation effect in the hemodynamic response.

The initial hemoglobin overshoot observed in human is not supported by the findings in small animals, for which the habituation effect in hemoglobin response appears for longer stimulus trains (8–20 s) (Martindale et al., 2005; Franceschini et al., 2008). This discrepancy may be due to differences between species (humans versus rats) or anticipation effects, which are suppressed by anesthesia in rats.

In addition to linearity/non-linearity of the responses, we performed a convolution analysis to test whether individual components of the current amplitude can predict the hemodynamic response time courses better than input stimulus. To do so, we used either a current amplitude component or the input stimulus as an input in the linear convolution model to predict the hemodynamic responses and determined statistically which factor shows better predictive power. With the stimulus tested, we found that there is no significant difference between the predictions of HbR using N20, P35 or P60 versus using the input stimulus. Instead, the sum of the peak amplitudes and the RMS value can predict HbR response consistently better than the input stimulus. Furthermore, subtracting N20 helps remove possible hyperpolarising DC shift of the membrane potential (Hellweg et al., 1977), both P35 − N20 and P60 − N20 achieve statistically better prediction of HbR than the stimulus.

The fact that the sum of three peaks and the RMS value of the current dipole amplitude are the best and N20 is the worst in predicting the hemodynamic responses is in agreement with our recent results where we investigated neurovascular coupling using EEG and DOI in rats (Franceschini et al., 2008). Similar to the human study, the animal study shows that the later peaks (N1 and P2) of the somatosensory evoked potentials are able to predict hemodynamic responses better than the first peak, P1. The EEG evoked potential deflection P1 in rats is equivalent to MEG N20 deflection in humans. This initial peak constitutes the primary cortical response which is generated by synaptic excitation of middle layers by thalamocortical inputs (Li et al., 1956; Mitzdorf, 1985; Di et al., 1990). In most invasive animal studies, this earliest peak is typically used to predict the hemodynamic response (Caesar et al., 2003; Jones et al., 2004; Sheth et al., 2004; Iadecola, 2004). Historically, the greater focus on P1 may be partially due to the fact that this peak is very stable and persists with deep anesthesia, while subsequent cortical activity is abolished with deep anesthesia. The later cortical activity is more spread both spatially and temporally, and it is less stable with changes in stimulus parameters (Cauller and Kulics, 1991). Our findings that the hemodynamic responses correlate better with later peaks than with the primary cortical response suggest that the cortical hemodynamic response is largely controlled by synaptic activity related to intracortical processing rather than direct thalamic input. This is probably because the majority of the synapses are triggered with column inputs from other neurons in the cortex, and the hemodynamic response is mostly sensitive to the level of synaptic activity (Mathiesen et al., 1998). Thus, our results indicate that later components need to be considered to correctly understand the neurovascular coupling.

Despite the better predictions of hemoglobin response by the sum of two or three peaks and the RMS value of the current dipole amplitude, the habituation effects exhibited in the hemodynamic response are not totally accounted for with these parameters. This suggests that the hemodynamic response cannot be described by a single component, but by a weighted combination of multiple components or some later components in the neural signals. A more complex experimental design which provides more variability in the neuronal responses will be necessary to obtain sufficient power in the statistical analysis of the different neurovascular coupling models.


We thank Dr. Solomon Diamond and Dr. Ted Huppert for helping with probe design and Deirdre Foxe and Natsuko Mori for helping with MEG measurements. This research is supported by the US National Institutes of Health (NIH) grants R01-EB001954, P41-RR14075 and R01-EB006385. Wanmei Ou is partially supported by the NSF graduate fellowship and the PHS training grant DA022759-03.

Appendix A. Supplementary data

Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.neuroimage.2009.03.008.


1For subjects 2 and 5, the window is set as 30–160 ms after each stimulus onset in a train; for the rest of the subjects, the window is set as 30–140 ms.

2The time windows of yHbO(t) and yHbR(t) for the linear convolution analysis are separately selected based on the zero crossing of the corresponding responses with respect to the baseline for each condition and each subject.


  • Allison T. Scalp and cortical recordings of initial somatosensory cortex activity to median nerve stimulation in man. Ann. NY. Acad. Sci. 1992;388:671–678. [PubMed]
  • Baillet S, Mosher JC, Leahy RM. Electromagnetic brain mapping. IEEE Signal Process. Mag. 2001;18:14–30.
  • Bandettini PA, Ungerleider LG. From neuron to BOLD: new connections. Nat. Neurosci. 2001;4:864–866. [PubMed]
  • Bandettini PA, Birn RM, Saad ZS, Kelley D. Dynamic nonlinearities in BOLD contrast: neuronal or hemodynamic? Int. Congr. Ser. 2002;1235:73–85.
  • Birn RM, Saad ZS, Bandettini PA. Spatial heterogeneity of the nonlinear dynamics in the fMRI BOLD response. NeuroImage. 2001;14:817–826. [PubMed]
  • Birn RM, Bandettini PA. The effect of stimulus duty cycle and “off” duration on BOLD response linearity. NeuroImage. 2005;27:70–82. [PubMed]
  • Boynton GM, Engel SA, Glover GH, Heeger DJ. Linear systems analysis of functional magnetic resonance imaging in human V1. J. NeuroSci. 1996;16:4207–4221. [PubMed]
  • Brenner D, Lipton J, Kaufman L, Williamson SJ. Somatically evoked magnetic fields of the human brain. Science. 1978;199:81–83. [PubMed]
  • Butti M, Pastori A, Merzagora A, Zucca C, Bianchi A, Reni G, Cerutti S. Multimodal analysis of a sustained attention protocol: Continuous Performance Test assessed with Near Infrared Spectroscopy and EEG. Conf. Proc. IEEE. Eng. Med. Biol. Soc. 2006;1:1040–1043. [PubMed]
  • Cauller LJ, Kulics AT. The neural basis of the behaviorally relevant N1 component of the somatosensory-evoked potential in SI cortex of awake monkeys: evidence that backward cortical projections signal conscious touch sensation. Exp. Brain Res. 1991;84:607–619. [PubMed]
  • Caesar K, Thomsen K, Lauritzen M. Dissociation of spikes, synaptic activity, and activity-dependent increments in rat cerebellar blood flow by tonic synaptic inhibition. Proc. Natl. Acad. Sci. U. S. A. 2003;100:16000–16005. [PubMed]
  • Dale AM, Sereno M. Improved localization of cortical activity by combining EEG and MEG with MRI cortical surface reconstruction: a linear approach. J. Cogn. Neurosci. 1993;5:162–176. [PubMed]
  • Dale AM. Optimal experimental design for event-related fMRI. Hum. Brain Mapp. 1999;8:109–114. [PubMed]
  • Dale AM, Fischl B, Sereno MI. Cortical surface-based analysis: I. segmentation and surface reconstruction. NeuroImage. 1999;9:179–194. [PubMed]
  • Dale AM, Liu AK, Fischl B, Buckner RL, Belliveau JW, Lewine JD, Halgren E. Dynamic statistical parametric mapping: combining fMRI and MEG for high-resolution imaging of cortical activity. Neuron. 2000;26:55–67. [PubMed]
  • Devor A, Dunn AK, Andermann ML, Ulbert I, Boas DA, Dale AM. Coupling of total hemoglobin concentration, oxygenation, and neural activity in rat somato-sensory cortex. Neuron. 2003;39:353–359. [PubMed]
  • Di S, Baumgartner C, Barth DS. Laminar analysis of extracellular field potentials in rat vibrissa/barrel cortex. J. Neurophysiol. 1990;63:832–840. [PubMed]
  • Fischl B, Sereno M, Dale AM. Cortical surface-based analysis: II. inflation, flattening, and a surface-based coordinate system. NeuroImage. 1999;9:195–207. [PubMed]
  • Franceschini MA, Fantini S, Thompson JH, Culver JP, Boas DA. Hemodynamic evoked response of the sensorimotor cortex measured noninvasively with near-infrared optical imaging. Psychophysiology. 2003;40:548–560. [PMC free article] [PubMed]
  • Franceschini MA, Nissilä I, Wu W, Diamond SG, Bonmassar G, Boas DA. Coupling between somatosensory evoked potentials and hemodynamic response in the rat. NeuroImage. 2008;41:189–203. [PMC free article] [PubMed]
  • Goff WR, Allison T, Vaughan HG. The functional neuroanatomy of event related potentials. In: Callaway E, Tueting P, Koslow SH, editors. Event-related Brain Potentials in Man. New York: Academic Press; 1978. pp. 1–79.
  • Hämäläinen MS, Ilmoniemi R. Interpreting Measured Magnetic Fields of the Brain: Estimates of Current Distributions. Helsinki University Technical Report TKK-F-A559. 1984
  • Hämäläinen MS, Sarvas J. Realistic conductivity geometry model of the human head for interpretation of neuromagnetic data. IEEE Trans. Biomed. Engr. 1989;36:165–171. [PubMed]
  • Hämäläinen MS. NMR Center. Harvard University: Mass General Hospital; 2005. MNE Software User's Guide.
  • Hari R, Reinikainen K, Kaukoranta E, Hämäläinen M, Ilmoniemi R, Penttinen A, Salminen J, Teszner D. Somatosensory evoked cerebral magnetic fields from SI and SII in man. Electroencephalogr. Clin. Neurophysiol. 1984;57:254–263. [PubMed]
  • Hari R, Forss N. Magnetoencephalography in the study of human somatosensory cortical processing. Philos. Trans. R. Soc. Lond. B. Biol. Sci. 1999;354:1145–1154. [PMC free article] [PubMed]
  • Hellweg F, Schultz W, Creutzfeldt O. Extracellular and intracellular recordings from cat's cortical whisker projection area: thalamocortical response transformation. J. Neurophysiol. 1977;40:463–479. [PubMed]
  • Herrmann MJ, Huter T, Plichta MM, Ehlis AC, Alpers GW, Mühlberger A, Fallgatter AJ. Enhancement of activity of the primary visual cortex during processing of emotional stimuli as measured with event-related functional near-infrared spectroscopy and event-related potentials. Human. Brain Mapp. 2008;29:28–35. [PubMed]
  • Iadecola C. Neurovascular regulation in the normal brain and in Alzheimer's disease. Nat. Rev. Neurosci. 2004;5:347–360. [PubMed]
  • Jones M, Berwick J, Johnston D, Mayhew J. Concurrent optical imaging spectroscopy and laser-Doppler flowmetry: the relationship between blood flow, oxygenation, and volume in rodent barrel cortex. Neuroimage. 2001;13:1002–1015. [PubMed]
  • Jones M, Hewson-Stoate N, Martindale J, Redgrave P, Mayhew J. Nonlinear coupling of neural activity and CBF in rodent cortex. Neuroimage. 2004;22:956–965. [PubMed]
  • Kaufman L, Okada Y, Brenner D, Williamson SJ. On the relation between somatic evoked potentials and fields. Int. J. Neurosci. 1981;15:223–239. [PubMed]
  • Koch SP, Koendgen S, Bourayou R, Steinbrink J, Obrig H. Individual alpha-frequency correlates with amplitude of visual evoked potential and hemodynamic response. NeuroImage. 2008;41:233–242. [PubMed]
  • Li CL, Cullen C, Jasper HH. Laminar microelectrode studies of specific somatosensory cortical potentials. J. Neurophysiol. 1956;19:111–130. [PubMed]
  • Logothetis NK, Wandell BA. Interpreting the BOLD signal. Annu. Rev. Physiol. 2004;66:735–769. [PubMed]
  • Mackert BM, Wubbeler G, Leistner S, Uludag K, Obrig H, Villringer A, Trahms L, Curio G. Neurovascular coupling analyzed non-invasively in the human brain. Neuroreport. 2004;15:63–66. [PubMed]
  • Mackert BM, Leistner S, Sander T, Liebert A, Wabnitz H, Burghoff M, Trahms L, Macdonald R, Curio G. Dynamics of cortical neurovascular coupling analyzed by simultaneous DC-magnetoencephalography and time-resolved near-infrared spectroscopy. Neuroimage. 2008;39:979–986. [PubMed]
  • Mantini D, Perrucci MG, Cugini S, Ferretti A, Romani GL, Del Gratta C. Complete artifact removal for EEG recorded during continuous fMRI using independent component analysis. NeuroImage. 2007;34:598–607. [PubMed]
  • Martindale J, Berwick J, Martin C, Kong Y, Zheng Y, Mayhew J. Long duration stimuli and nonlinearities in the neural-haemodynamic coupling. J. Cereb. Blood Flow Metab. 2005;25:651–661. [PubMed]
  • Mathiesen C, Caesar K, Akgören N, Lauritzen M. Modification of activity-dependent increases of cerebral blood flow by excitatory synaptic activity and spikes in rat cerebellar cortex. J. Physiol. 1998;512:555–566. [PubMed]
  • Mauguiere F, Merlet I, Forss N, Vanni S, Jousmäki V, Adeleine P, Hari R. Activation of a distributed somatosensory cortical network in the human brain. A dipole modeling study of magnetic fields evoked by median nerve stimulation. Part I: Location and activation timing of SEF sources. Electroencephalogr. Clin. Neurophysiol. 1997;104:281–289. [PubMed]
  • Mitzdorf U. Current source-density method and application in cat cerebral cortex: investigation of evoked potentials and EEG phenomena. Physiol. Rev. 1985;65:37–100. [PubMed]
  • Mosher JC, Lewis PS, Leahy RM. Multiple dipole modeling and localization from spatio-temporal MEG data. IEEE Trans. Biomed. Engr. 1992;39:541–557. [PubMed]
  • Negishi M, Abildgaard M, Nixon T, Constable RT. Removal of time-varying gradient artifacts from EEG data acquired during continuous fMRI. Clin. Neurophysiol. 2004;115:2181–2192. [PubMed]
  • Nelder JA, Mead R. A simplex method for function minimization. Comput. J. 1965;7:308–313.
  • Ngai AC, Jolley MA, D'Ambrosio R, Meno JR, Winn HR. Frequency-dependent changes in cerebral blood flow and evoked potentials during somatosensory stimulation in the rat. Brain. Res. 1999;837:221–228. [PubMed]
  • Norup Nielsenq A, Lauritzen M. Coupling and uncoupling of activity-dependent increases of neuronal activity and blood flow in rat somatosensory cortex. J. Physiol. 2001;533:773–785. [PubMed]
  • Obrig H, Neufang M, Wenzel R, Kohl M, Steinbrink J, Einhäupl K, Villringer A. Spontaneous low frequency oscillations of cerebral hemodynamics and metabolism in human adults. Neuroimage. 2000;12:623–639. [PubMed]
  • Obrig H, Israel H, Kohl-Bareis M, Uludag K, Wenzel R, Müller B, Arnold G, Villringer A. Habituation of the visually evoked potential and its vascular response: implications for neurovascular coupling in the healthy adult. NeuroImage. 2002;17:1–18. [PubMed]
  • Oostendorp TF, van Oosterom A. Source parameter estimation in inhomogeneous volume conductors of arbitrary shape. IEEE Trans. Biomed. Engr. 1989;36:382–391. [PubMed]
  • Riera JJ, Jimenez JC, Wan X, Kawashima R, Ozaki T. Nonlinear local electrovascular coupling II: from data to neuronal masses. Human. Brain Mapp. 2007;28:225–254. [PubMed]
  • Ritter P, Villringer A. Simultaneous EEG-fMRI. Neurosci. Biobehav. Rev. 2006;30:823–838. [PubMed]
  • Rovati L, Salvatori G, Bulf L, Fonda S. Optical and electrical recording of neural activity evoked by graded contrast visual stimulus. Biomed. Eng. Online. 2007;6:28. [PMC free article] [PubMed]
  • Royl G, Leithner C, Sellien H, Müller JP, Megow D, Offenhauser N, Steinbrink J, Kohl-Bareis M, Dirnagl U, Lindauer U. Functional imaging with Laser Speckle Contrast Analysis: vascular compartment analysis and correlation with Laser Doppler Flowmetry and somatosensory evoked potentials. Brain Res. 2006;1121:95–103. [PubMed]
  • Sander TH, Liebert A, Mackert BM, Wabnitz H, Leistnerz S, Curio G, Burghoff M, Macdonald R, Trahms L. DC-magnetoencephalography and time-resolved near-infrared spectroscopy combined to study neuronal and vascular brain responses. Physiol. Meas. 2007;28:651–664. [PubMed]
  • Sato H, Kiguchi M, Kawaguchi F, Maki A. Practicality of wavelength selection to improve signal-to-noise ratio in near-infrared spectroscopy. NeuroImage. 2004;21:1554–1562. [PubMed]
  • Schmid MC, Oeltermann A, Juchem C, Logothetis NK, Smirnakis SM. Simultaneous EEG and fMRI in the macaque monkey at 4.7 Tesla. Magn. Reson. Imaging. 2006;24:335–342. [PubMed]
  • Sheth SA, Nemoto M, Guiou M, Walker M, Pouratian N, Toga AW. Evaluation of coupling between optical intrinsic signals and neuronal activity in rat somatosensory cortex. NeuroImage. 2003;19:884–894. [PubMed]
  • Sheth SA, Nemoto M, Guiou M, Walker M, Pouratian N, Toga AW. Linear and nonlinear relationships between neuronal activity, oxygen metabolism, and hemodynamic responses. Neuron. 2004;42:347–355. [PubMed]
  • Source modeling software (xfit) Finland: Elekta-Neuromag Oy, Helsinki;
  • Strangman G, Franceschini MA, Boas DA. Factors affecting the accuracy of near-infrared spectroscopy concentration calculations for focal changes in oxygenation parameters. Neuroimage. 2003;18:865–879. [PubMed]
  • Tiihonen J, Hari R, Hämäläinen M. Early deflections of cerebral magnetic responses to median nerve stimulation. Electroencephalogr. Clin. Neurophysiol. 1989;74:290–296. [PubMed]
  • Uludag K, Steinbrink J, Villringer A, Obrig H. Separability and cross talk: optimizing dual wavelength combinations for near-infrared spectroscopy of the adult head. NeuroImage. 2004;22:583–589. [PubMed]
  • Ureshi M, Matsuura T, Kanno I. Stimulus frequency dependence of the linear relationship between local cerebral blood flow and field potential evoked by activation of rat somatosensory cortex. Neurosci. Res. 2004;48:147–153. [PubMed]
  • Vasios CE, Angelone LM, Purdon PL, Ahveninen J, Belliveau JW, Bonmassar G. EEG/(f)MRI measurements at 7 Tesla using a new EEG cap (“InkCap”) NeuroImage. 2006;33:1082–1092. [PubMed]
  • Vazquez AL, Noll DC. Nonlinear aspects of the BOLD response in functional MRI. Neuroimage. 1998;7:108–118. [PubMed]
  • Wikström H, Huttunen J, Korvenoja A, Virtanen J, Salonen O, Aronen H, Ilmoniemi RJ. Effects of interstimulus interval on somatosensory evoked magnetic fields (SEFs): a hypothesis concerning SEF generation at the primary sensorimotor cortex. Electroencephalogr. Clin. Neurophysiol. 1996;100:479–487. [PubMed]
  • Wray S, Cope M, Delpy DT. Characteristics of the near infrared absorption spectra of cytochrome aa3 and hemoglobin for the noninvasive monitoring of cerebral oxygenation. Biochim. Biophys. Acta. 1988;933:184–192. [PubMed]
  • Yamashita Y, Maki A, Kiozumi H. Wavelength dependence of the precision of noninvasive optical measurement of oxy-, deoxy-, and total-hemoglobin concentration. Med. Phys. 2001;28:1108–1114. [PubMed]
  • Zhang Y, Brooks DH, Boas DA. A haemodynamic response function model in spatio-temporal diffuse optical tomography. Phys. Med. Biol. 2005;50:4625–4644. [PubMed]