Modeling Neurovascular Coupling from Clustered Parameter Sets for Multimodal EEG-NIRS

Talukdar, M. Tanveer; Frost, H. Robert; Diamond, Solomon G.

doi:https://doi.org/10.1155/2015/830849

Computational and Mathematical Methods in Medicine

On this page

Abstract Introduction Methods Results and Discussion Conclusion Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2015 | Article ID 830849 | https://doi.org/10.1155/2015/830849

Modeling Neurovascular Coupling from Clustered Parameter Sets for Multimodal EEG-NIRS

M. Tanveer Talukdar,¹H. Robert Frost,²and Solomon G. Diamond¹

Academic Editor: Phaneendra K. Yalavarthy

Received24 Nov 2014

Revised03 Feb 2015

Accepted04 Feb 2015

Published19 May 2015

Abstract

Despite significant improvements in neuroimaging technologies and analysis methods, the fundamental relationship between local changes in cerebral hemodynamics and the underlying neural activity remains largely unknown. In this study, a data driven approach is proposed for modeling this neurovascular coupling relationship from simultaneously acquired electroencephalographic (EEG) and near-infrared spectroscopic (NIRS) data. The approach uses gamma transfer functions to map EEG spectral envelopes that reflect time-varying power variations in neural rhythms to hemodynamics measured with NIRS during median nerve stimulation. The approach is evaluated first with simulated EEG-NIRS data and then by applying the method to experimental EEG-NIRS data measured from 3 human subjects. Results from the experimental data indicate that the neurovascular coupling relationship can be modeled using multiple sets of gamma transfer functions. By applying cluster analysis, statistically significant parameter sets were found to predict NIRS hemodynamics from EEG spectral envelopes. All subjects were found to have significant clustered parameters () for EEG-NIRS data fitted using gamma transfer functions. These results suggest that the use of gamma transfer functions followed by cluster analysis of the resulting parameter sets may provide insights into neurovascular coupling in human neuroimaging data.

1. Introduction

Neural activity is a complex biophysical process that involves electrochemical and vascular interaction at the cellular level. Modern neuroimaging technologies such as functional magnetic resonance imaging (fMRI), near infrared spectroscopy (NIRS), and diffuse optical imaging (DOI) typically detect changes in the local blood flow associated with the vascular response to provide an indirect measure of neural activity. A major limitation with these techniques is their inability to map precisely the source of the underlying neural activity [1, 2]. Consequently, understanding neurovascular coupling or the relationship linking the underlying neural activity to the local changes in cerebral hemodynamics remains a vital area of research.

To date, noninvasive multimodal neuroimaging techniques with electroencephalography (EEG) and fMRI have been widely used to investigate the relationship between neural activity and cerebral hemodynamics [3–5]. These studies attempt to show quantitative aspects such as the degree of correlation between electrical scalp potentials recorded using EEG and the hemodynamic response measured from fMRI after separately preprocessing the acquired data. Only a few studies have investigated the direct link between neural activity and the hemodynamic response by estimating transfer functions which map neural rhythms recorded using EEG to fMRI blood-oxygen level dependent (BOLD) signals reflecting the hemodynamic response. Despite being informative, there are a number of methodological and technical difficulties associated with multimodal EEG-fMRI techniques. EEG data recorded during fMRI acquisitions can be contaminated by gradient artifacts induced by the changing magnetic field gradients used for spatial encoding in MRI [6]. Electromagnetic interference could decrease the signal quality in both modalities and therefore obscure features of interest in the datasets being used to identify the neurovascular coupling relationship. Moreover, EEG-fMRI instrumentation can be fairly complex to set up and is not readily accessible for conducting neuroimaging research in smaller research facilities.

Multimodal neuroimaging with EEG and NIRS is increasingly gaining popularity due to the portability of EEG-NIRS instrumentation. EEG-NIRS systems can provide measurements with high signal-to-noise ratio (SNR) since there is no electrooptical interference. More importantly, NIRS can measure the hemodynamic changes at the capillary level while fMRI detects BOLD signal that is sensitive to changes in venous outflow [7]. This makes NIRS potentially advantageous for measuring hemodynamics associated with neural activity as there is evidence which suggests that neuronal function is supported by brain capillary oxygenation [8]. A number of studies have used EEG-NIRS systems to investigate neurovascular coupling from the measurement of the electrical scalp potential and the vascular response in stimuli evoked neural responses [9–11]. Although these prior works show that neural correlates of EEG signals covary with NIRS hemodynamics, there is still not a clear understanding of how the signals measured with EEG are predictive of NIRS. Newer methods need to be established based on models that relate neural activity to the hemodynamic response.

In this work we introduce a neurovascular-coupling modeling approach based on fitting gamma transfer functions and then clustering the parameters into sets that most effectively map the electrical scalp potentials recorded using EEG to the cerebrovascular response measured from NIRS hemodynamics. Gamma transfer functions have been used to represent the hemodynamic response function (HRF) in fMRI analysis methods [12–14]. Typically in fMRI data analysis, the BOLD signal is regressed on a general linear model (GLM) constructed by convolving the HRF with a boxcar function representing the temporal structure of the experimental paradigm [15]. In recent years, NIRS studies have also adopted the GLM framework to map brain activations [16]. Unfortunately, a boxcar function convolved with a canonical HRF may not accurately represent the cerebrovascular response since there is variability in the neural response between brain regions and different tasks [17, 18]. In addition, the shape of the HRF will depend on the type of imaging modality used [19]. NIRS hemodynamics, which have higher temporal resolution than fMRI BOLD signals and different hemodynamic measures, will require a different parameter set for the HRF compared to parameters used in the BOLD fMRI. With respect to the GLM modeling framework, our neurovascular-coupling modeling approach takes advantage of the measured EEG signal linked directly to neural rhythms in order to predict the NIRS hemodynamic response. Moreover, we vary the parameters of the gamma transfer function so that our model predictions have improved fit. Finally, using our approach we want to identify reliable gamma transfer function parameters that can predict NIRS HbO hemodynamic response from EEG spectral envelopes within an individual.

2. Methods

We conducted an Institutional Review Board (IRB) approved experiment that involved unilateral stimulation of the median nerve. Typically, stimulation of the median nerve on the hand induces mu rhythm desynchronization in the 8–13 Hz range in the contralateral somatosensory area of the brain. The cortical response can be detected as fluctuations in the EEG mu band and is also known to produce positive changes in HbO calculated from NIRS optical measurements [20, 21]. Our experimental paradigm is designed in order to facilitate detection of signal components that correlate with the median nerve stimulus. In this way we have selected for signals that are likely to have originated in the same brain region, which makes it practical to study neurovascular coupling.

In our study we use an EEG-NIRS head probe that was designed to record neural activity from optimally placed EEG electrodes and NIRS optodes, which cover the whole head [22]. The positioning of the EEG sensors is configured so that they are equidistant and close in proximity to a NIRS optode. This ensures to some degree that the measured signals from an EEG electrode and a NIRS optode have a common source in the cerebral cortex. The EEG-NIRS instrumentation consists of 16 optical sources and 8 detectors interleaved with EEG sensors, which are arranged according to the International 10–20 system. The NIRS system emits continuous wave (CW) near infrared light at four wavelengths (690 nm, 785 nm, 808 nm, and 830 nm). The source-detector spacing is 33 mm, which allows for probing neural activation with high sensitivity [23]. The EEG system measures the electrical scalp potentials from 30 electrodes. The data acquisition rates for the EEG and NIRS systems are 2048 Hz and 25 Hz, respectively.

2.1. Data Preprocessing

EEG data is contaminated by ocular, muscular, cardiac, and other physiological artifacts. For the sampled EEG data, epochs containing aberrant waveforms were manually discarded by visual inspection using the EEGLAB TOOLBOX 10.2.2.4b (Swartz Center for Computational Neurosciences, La Jolla, CA; http://www.sccn.ucsd.edu/eeglab/). We then applied independent component analysis (ICA) using the Infomax ICA algorithm [24], as implemented in EEGLAB to remove ocular movement and blink artifacts from the EEG data [25–27]. Empirical mode decomposition (EMD) was next applied to the motion corrected and noise artifact reduced EEG data for extracting neural rhythms linked to mu rhythm desynchronization. EMD decomposes a signal into a set of band limited functions called intrinsic oscillatory modes or integral mode functions (IMF) [28]. Each IMF represents oscillations in a narrow frequency band and may reveal particular features of interest in the given signal. The EMD method can offer better time-frequency localization compared to traditional methods like short time Fourier transform (STFT) and wavelet decomposition techniques [29, 30]. STFT has fixed time-frequency resolution due to the fact that it gives a global frequency distribution in the processing time window. Wavelet analysis, on the other hand, can provide variable time-frequency resolution compared to STFT method. However, there is a trade off between the temporal and spatial scales depending on choice of the wavelet basis function. In EMD, one does not require any basis functions and the signals do not need to be stationary, as is assumed in STFT and wavelet analysis methods. As a result EMD has gained considerable popularity for the analysis of EEG data, which exhibit a high degree of nonlinearity and oscillatory rhythms in subband frequencies [31, 32].

NIRS optical signal can be contaminated by motion artifacts and sources of physiological noise such as cardiac signals. Motion artifacts were rejected by applying Chauvenet’s criterion [33]. A deviation ratio (DR) at each time point was calculated by dividing the NIRS signal deviations by the standard deviation of the signal deviations. Signal deviations were calculated as the difference between the raw data and its smoothed version obtained by applying a moving average filter with a span of 30 points (1.2 s duration). Data with a DR greater than the standard Chauvenet’s criterion threshold were eliminated and the discontinuous data segments then spliced together. The resulting NIRS signal was then used to calculate changes in hemodynamics by solving equations of the modified Beer-Lambert law (MBLL), which relates the attenuation of light to relative changes in concentration of oxy-hemoglobin (HbO) and deoxy-hemoglobin (HbR) [34]. HbO and HbR were then detrended and low pass filtered using a zero-phase 3rd order Butterworth filter with a cutoff frequency of 0.5 Hz to remove cardiac oscillations, which contaminate NIRS signal [35].

Figure 1 shows a flowchart of the preprocessing stream for the simultaneously acquired EEG-NIRS data during median nerve stimulation. EEG and NIRS data were sampled at 2048 Hz and 25 Hz, respectively, using the prototype EEG-NIRS head probe. Both EEG and NIRS data were treated for artifact rejection and were matched temporally across sampled data points. The EEG recordings were bandpass filtered (5–14 Hz) using a zero-phase 4th order Butterworth filter in order to contain the mu rhythm band (8–13 Hz). Hilbert transform [36] was then applied on the first IMF signal component (IMF1) to generate EEG spectral envelopes, which was later downsampled to 25 Hz. The EEG spectral envelope corresponding to IMF1 was chosen because it is modulated by frequency components in the mu band and is temporally correlated to the NIRS hemodynamic response due to the median nerve stimuli. The EEG spectral envelope and NIRS hemodynamics also have similar timing and bandwidth in the frequency domain.

2.2. Neurovascular-Coupling Model

In our model of neurovascular coupling, we assume that gamma transfer functions can map EEG spectral envelopes containing neural rhythms to NIRS HbO reflecting the hemodynamic response linked to the underlying neural activity. This mapping can be expressed bywhere is an EEG spectral envelope, which is convolved with a candidate gamma transfer function to predict the hemodynamic response . The variables and are the gain and offset, respectively. Gamma transfer functions can be expressed in exponential form with parameters , , and [37]: where determines the rise time to peak amplitude, is an integer which governs the shape of the function, and is a pure delay. In our method, we recover the parameters , , and of the gamma transfer function that produces the best fit between the observed and predicted NIRS hemodynamics. Identification of the gamma transfer function and its parameters are carried out in two steps. In the first step, a brute force technique is applied to calculate values and in (1) by least squares inversion for each gamma transfer function constructed from a set of predefined parameters , , and . These predefined parameter values ( = 0.1–0.6 s, = 1–4, and = 0.5–3 s) were empirically chosen based on the observation of NIRS HbO response time to median nerve stimulation. Usually the delay in HbO response time can vary between 2 and 3 s, while the time to peak can be fairly quick 0.1–0.5 s [38]. From the solution space, we then identify , , , , and , which minimize the residual sum squared error (SSE). In the second step, the identified gamma function parameters (, , and ) including the gain and offset are jointly optimized using a simplex search method [39]. This is performed using Maltab’s “fminsearch” routine in which we specify SSE as the objective function.

2.3. Estimating Gamma Transfer Functions for Simulated EEG-NIRS Data

We simulated EEG data reflecting mu rhythm desynchronization associated with neural responses in the somatosensory cortex [40]. The artificially produced EEG data was in the form of an amplitude modulated sinusoidal signal having a carrier frequency of 9 Hz. A rectangular pulse train with eight pulses having a width and period of 20 s and 40 s, respectively, was used to modulate the sinusoidal signal. This pulse train was first convolved with a Kaiser window having a size of 250 data points. This smoothed the edges of the pulse train and introduced a time shift of approximately 5 s in the pulse onset times. The sign of the time shifted pulse train was then reversed to represent the spectral envelope for mu rhythm desynchronization events that diminish in power during the stimulus. This was then used to generate the simulated EEG signal by standard amplitude modulation operation given by where is the modulation wave, is the carrier wave, and is the modulated signal. The EEG signal thus contained eight 20 s mu rhythm desynchronization events with a latency of 5 s at each onset time.

NIRS HbO was generated from the forward model of (1). Gamma transfer functions were evaluated for five different combinations of , , and that ranged in values from 0.3–0.6 s, 2–4, and 2-3 s, respectively. The gain and offset were set to 1 and 0. Each of the 5 gamma transfer functions was next convolved with the amplitude signal of the first IMF (IMF1) obtained from decomposition of the simulated EEG by the EMD technique. IMF1 basically had a spectral peak around the carrier frequency of 9 Hz and its amplitude signal was calculated from the Hilbert transform [36]. Both the lMF signal and NIRS HbO were normalized in their amplitudes by dividing by their respective standard deviations. We then titrated both the EEG IMF1 signal and NIRS HbO with varying degrees of white Gaussian noise (WGN) at each epoch. In case of HbO, we also added a sinusoidal signal component of 0.5 Hz frequency modulating the WGN. The resulting HbO were representative of more realistic NIRS measurements, which contain physiological noise such as cardiac oscillations that are in the 0.4–2.0 Hz range [41]. The amplitude ratio (AR) of the WGN was selected to be in the range from dB (no noise) to 16 dB.

After generating the 5 EEG-NIRS datasets, we estimated the known gamma transfer functions that were used to model the NIRS data and recover the corresponding parameters , , and . We applied the brute force and the simplex search method to calculate the gamma transfer function parameters for each EEG-NIRS epoch data segments that spanned the pulse start times to 15 s past their end times. The parameters , , and identified at each epoch data segments were then used in the forward model of (1) to predict the simulated HbO. The Pearson’s correlation between the simulated and predicted HbO was then evaluated at each amplitude ratio of the WGN.

2.4. Estimating Gamma Transfer Functions for Experimental EEG-NIRS Data

Gamma transfer functions that provide insight into the neurovascular coupling relationship in humans were estimated from EEG-NIRS data measured from three right handed male subjects (mean age of 31) in an IRB approved study. The subjects participating in the study had no history of neurological or psychiatric disorders and were in general good health. The left and right median nerves were stimulated using an electrical pulse generator that operated at 10 Hz. Pulses were transmitted during a 15-second block interval followed by 30 second rest interval during each experimental run that lasted a total of 6 minutes. Five recording sessions were carried out, two for each of the left and right median nerve stimulation runs and and one for which the subject was at rest and no stimulation was applied.

IMF1 spectral envelopes were derived for each of the 30 EEG channels. We then selected sets of IMF1 spectral envelopes corresponding to EEG channels located in the left and right hemispheres based on how strongly they correlated with the electrical pulse sequence. Similarly, NIRS HbO signals were selected within each hemisphere if they showed changes in concentrations above 0.5 μM. For each EEG IMF1 spectral envelope and NIRS HbO paired dataset within a hemisphere, we extracted epochs spanning the stimulus interval of 15 s including 5 s of pre- and poststimulus period. The epoch data were later smoothed with a moving average window having a time span of 4 s and normalized by dividing by their standard deviation. Each of the epoch datasets consisting of normalized IMF1 spectral envelope and NIRS HbO finally served as the neurovascular-coupling model input and output, respectively. Gamma transfer functions were then estimated for each epoch data segment using the two-step method described in Section 2.2.

We next investigated whether the estimated gamma transfer functions could be clustered into groups which are statistically significant in their predictions of the NIRS HbO hemodynamics. By clustering the estimated gamma transfer functions, we hope to reduce the effect of overfitting in our neurovascular-coupling model. In addition, the clustering method is intended to identify reliable gamma transfer functions which link the EEG input to the NIRS hemodynamic output at the individual level. For each of the estimated gamma transfer functions, the parameters , , , gain , and offset were used to form a 5-dimensional feature vector or data point. Since the features were scaled differently, we standardized the feature space so that they have zero mean unit variance. Hierarchical clustering was then applied on the feature space derived from all the standardized parameter feature vectors comprised of , , , , and . We chose this clustering method due to its relative ease of implementation and also because it is widely used as a clustering tool. In hierarchical clustering, data is grouped by linking feature vectors in a binary tree called a dendrogram [42]. There are a number of linkage methods to merge the feature vectors into clustered parameter sets. We used Ward’s linkage method as it has been shown to have better clustering performance [43]. In this method the within-cluster variance is minimized over all partitions obtainable by merging two clusters from the previous generation.

By specifying threshold values for cutting off the dendrogram at specific depths, we clustered the feature vectors into cluster divisions , where . For each clustered parameter set, we calculated the Pearson’s correlation coefficient between the predicted NIRS HbO obtained using the forward model (Equation (1)) and the measured HbO. The correlations were Fisher transformed to give values that are approximately standard normal in their distribution. The mean value, , was then computed for each of the clustered parameter sets. The value was next tested for statistical significance by performing a two-sided -test under the null hypothesis that = 0, which corresponds to a correlation coefficient of zero. This was done by first computing the score, , from the relationship: where is the number of feature vectors in the cluster. This value was then compared against the standard score at the 5% significance level (). It is important to note that the values generated via this two-sided -test are only approximate, are likely biased towards more extreme values, and should not be interpreted as accurate measures of statistical significance. This inaccuracy is due to both the fact that the underlying EEG-NIRS data does not follow independent observations from a bivariate normal distribution under the null hypothesis, giving the Fisher-transformed Pearson correlation coefficients a distribution different from , and the fact that -statistics computed for each gamma transfer function parameter vector within a given cluster are not independent, increasing the variance of the mean of the -statistics. Despite this inaccuracy, the computed values can still be successfully used to help to select the optimal split level for the dendrogram and to identify clusters of gamma transfer function parameter vectors that closely model the empirical EEG-NIRS relationship and merit more detailed analysis and interpretation.

Due to the fact that we are performing significance tests on values in cluster divisions, the values were adjusted for multiple hypothesis correction (MHC) using the Bonferroni correction [44]. For each cluster division 2 to 12, we also recorded the number of significant clusters by counting values that exceeded . From this count we assessed the optimal that the data can be divided into using hierarchical clustering. Our approach to MHC was to control the family-wise error rate (FWER) using the Bonferroni method for the family of hypothesis tests associated with the two-sided -tests performed for each of the clusters generated by one split of the dendrogram. Although MHC for this problem could have been used to control the false discovery rate (FDR) rather than the FWER or been applied to the larger family of hypothesis tests associated with all dendrogram splits, we choose the current approach to facilitate selection of the optimal dendrogram split level.

3. Results and Discussion

3.1. Simulated EEG-NIRS Data

Figure 2(a) shows an epoch of the simulated IMF1 spectral envelope (blue) and the pulse duration (green). The 5 gamma transfer functions that were used to simulate NIRS HbO are shown in Figure 2(b). The parameters of the gamma transfer function were chosen heuristically based on the notion that they reflect physiologically relevant NIRS hemodynamic response with respect to latency, shape, and peak time behavior.

(a)

(b)

Figures 3(a)–3(h) show simulated noise titrated HbO epoch data (red) obtained from convolution of the IMF1 spectral envelope in Figure 2(a) with a gamma transfer function having parameters s, , and s (see gamma function in blue, Figure 2(b)). For each of the 8 HbO epoch data, the estimated HbO are shown overlaid in black. Perfect fit is obtained when the noise AR is dB, which is practically zero noise level. At higher levels of noise, the estimated HbO has poorer fits with the simulated HbO. However, the estimated HbO in Figures 3(b)–3(d) seems to follow the dynamics of the noise-free simulated HbO in Figure 3(a). At noise AR of 9 dB and higher, the estimated HbO no longer exhibit the dynamics of noise-free simulated HbO. Similar trends are observed for HbO epoch data simulated using the other gamma transfer functions.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 4 shows the Pearson’s correlation values between the estimated HbO and the simulated HbO numbered 1 through 5 against noise AR in WGN. For the noise-free (AR = dB) simulated HbO numbered 1–5, = 1, indicating that the estimated and simulated HbO are in agreement with each other with a goodness of fit statistic . There is a sharp drop in value from 1 to 0.6 as AR increase from dB to 2.2 dB. At higher noise amplitude ratio from 2.2 dB to 11.5 dB, values drop almost linearly with a gentler slope. Between AR values of 12 dB to 16 dB, the values can be seen to level off. From this analysis we can note that reliable estimates of NIRS HbO, which is strongly correlated with the simulated HbO () can be obtained when the noise AR is below 4 dB.

Figure 5 shows parameters recovered at different AR values for the 5 sets of NIRS HbO generated from the known gamma transfer functions. We can note that, for AR = dB, the parameters recovered from the gamma transfer functions fits are identical to the original ones used to simulate HbO (see Figure 2(b)). However, all the recovered gamma transfer function parameters vary from the original ones for subsequent fits in which the AR values of the WGN increase from 2.2 dB to 16 dB. However, a number of the estimated gamma transfer function parameters remain stable in their value up to certain noise AR ratio. For example, the recovered parameters , , and for noise AR of 2.2 dB to 6.86 dB agree closely with the known gamma transfer function parameters with , , and (green dotted line). A similar observation can be made for the gamma transfer function with parameters , , and (orange dotted line) for which the parameter estimates of , , and agree closely with the known values as the noise AR increases from dB to 4.57 dB. It is also interesting to note that, in Figure 5(b), the shape parameter for 4 of the 5 gamma transfer functions seems to level off at around for noise AR greater than 6.86 dB. Despite the differences between the recovered and known values of the gamma transfer function parameter at increasing noise AR, our neurovascular-coupling model is still able to predict the dynamics of the simulated NIRS HbO containing moderate levels of WGN (see Figure 3).

(a)

(b)

(c)

3.2. Experimental EEG-NIRS Data

Figures 6(a)–6(d) show sample results from fitting experimental EEG-NIRS data from Subject 1. In each plot, the estimated NIRS HbO (purple) is overlaid on the calculated HbO (red). The EEG IMF1 spectral envelope is shown in blue and the pulse duration (15 s) in green. The EEG-NIRS data displayed in each plot have been normalized by their standard deviations. We can note that, after the onset of the stimulus at 0 s, the EEG IMF1 envelope decreases in amplitude while the NIRS HbO show a steady increase.

(a)

(b)

(c)

(d)

Figures 7(a)–7(k) shows cluster values computed for cluster division ranging from 2 to 12 in subject 1’s left hemisphere during right median nerve stimulation. The total number of feature vectors used in hierarchical clustering for this subject and experimental condition was 216. The statistical threshold for the two sided -test, which was adjusted for multiple comparison by Bonferroni correction, is indicated by the red line in each plot. The number of significant clusters with exceeding decreases with an increasing number of cluster divisions. For cluster divisions 11 and 12, there are no significant cluster values.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

Figure 8 shows the relationship between and , which were obtained from the plots shown in Figure 7. Since cluster divisions 2 to 5 have the highest number of significant clusters (), we take the optimal to be 5. This is based on the assumption that a higher number of clustered gamma transfer function parameters sets have greater flexibility for modeling the EEG-NIRS relationship compared to fewer clustered parameter sets.

Figure 9(a) shows the dendrogram obtained at optimal cluster division in subject 1 (Figure 8). The representative gamma transfer functions for the clustered parameters sets in the dendrogram are shown in Figure 9(b). Each of these functions was generated by taking the mean of the gamma transfer function parameters values corresponding to individual cluster groups. Their amplitudes were then scaled by dividing by the respective area under the curve. The blue and the green gamma transfer functions correspond to the cluster values that were statistically significant; that is, . The gamma transfer function in blue has a delay of about 2 s and peak time that is less than 0.5 s. In contrast, the gamma transfer function in green has a larger delay of approximately 3 s and a longer peak time of about 1 s.

(a)

(b)

Figure 10 shows subject 2 s statistically significant clusters for each cluster division . Hierarchical clustering was performed on feature vectors obtained by fitting EEG IMF1 envelopes to NIRS sequences in the left hemisphere during right median nerve stimulation. There were 96 feature vectors for this subject and experimental condition. Based on the graph of versus in Figure 10, the optimal cluster division is and the number of significant clusters is .

Figure 11(a) shows subject 2 s dendrogram for optimal cluster division identified in Figure 10. The representative gamma transfer functions corresponding to the three visible clusters in the dendrogram are shown in Figure 11(b). The gamma transfer function in red was identified as statistically significant. It has a delay of approximately 2.5 s and a peak time of about 1.5 s.

(a)

(b)

Table 1 shows the number of significant clusters identified for an optimal cluster division . The values shown correspond to the results from clustering recovered parameters of the gamma transfer function fit applied to EEG-NIRS data measured contralateral to the side of the median nerve stimulation.

4. Conclusion

Our mathematical model of neurovascular coupling is able to show the relationship between NIRS HbO and EEG spectral envelopes based on gamma transfer function fits with clustered parameter sets. In case of simulated EEG-NIRS data, the correlations between the predicted HbO and the simulated HbO are approximately 1 for zero WGN. The recovered gamma transfer function parameters (see Figure 5) are identical to the original ones in the absence of WGN (AR = dB). At higher noise amplitude ratios in WGN, the correlations drop almost linearly for all 5 NIRS datasets and we can still recover parameters of the gamma transfer function fits.

From the analysis of experimental EEG-NIRS data, we observed that multiple gamma transfer functions can significantly predict NIRS HbO hemodynamic response from EEG spectral envelopes representing mu rhythm desynchronization events. By clustering the recovered parameters of the gamma functions including the estimated gains and offsets, we found statistically significant gamma transfer functions that can predict NIRS HbO from EEG spectral envelopes. Those representative gamma transfer functions are characterized by different onset times and peak latencies in individual subjects (see Figures 9 and 11). These differences could be attributed to variability in the neural and hemodynamic response across brain regions including individual trials as suggested by previous studies [17, 19, 45].

This new method of analyzing the neurovascular coupling relationship has several potential advantages. It offers greater flexibility in modeling the relationship between EEG and NIRS data based on gamma transfer functions that can vary in their parameters. Using fixed transfer functions can introduce bias in modeling the cerebrovascular response and lead to inaccurate assessment of the onset times and latencies in neural activation [46, 47]. Also in fitting the gamma transfer functions, there are only five parameters to estimate. Hence, there is less potential for error in estimating them compared to models with a large number of parameters. Finally, we are using a data driven approach requiring minimal assumptions about the underlying biophysical process associated with neurovascular coupling. This makes the neurovascular coupling modeling problem more tractable and applicable in neuroimaging studies.

However, there may be some limitations in our approach to modeling the neurovascular coupling relationship. There is the potential for overfitting the data by allowing each gamma transfer function to have distinct parameters. We expect to overcome this problem by using the representative gamma transfer function derived from taking the average of the statistically significant parameter sets recovered from the model fitting of EEG-NIRS data measured from individual subjects. Due to the small number of human subjects, it is not possible to generalize about the nature of neurovascular coupling from the estimated parameters. However, the statistical testing was performed on an individual subject basis with a large number of feature vectors for each subject. This data was sufficient to allow us to develop this methodology. Nevertheless, acquiring EEG-NIRS data from a larger pool of subjects is necessary to validate neurovascular-coupling parameters that can be computed with this method.

Our method could be eventually extended to study the neurovascular coupling relationship in patients with neurodegenerative disorders. Several studies have reported that neurovascular coupling can be disrupted by brain disorders like Alzheimer’s disease (AD), stroke, and epilepsy [48]. Certain brain disorders may cause changes in the chemical mediators of neurovascular coupling, which could lead to abnormal patterns of vasodilation and prevent removal of harmful byproducts of molecular metabolism. Ischemic stroke in the brain for instance can limit adequate flow of oxygenated blood to active neurons leading to vasoparalysis and impairing brain function [49]. We hope that we can apply our analysis method to detect abnormal neurovascular coupling in clinical neuroimaging studies using multimodal EEG-NIRS systems. Based on the data collected from patients, it will be possible to identify clustered gamma function parameters that could be significantly different from normal healthy populations.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors thank both Dr. Paolo Giacometti and Dr. Katherine Purdue for help with the experimental setup and data acquisition using the EEG-NIRS head probe. Funding for this study has been provided by the Neukom Institute and the Institute for Quantitative Biomedical Scienes (IQBS) at Dartmouth.

References

B. Horwitz, K. J. Friston, and J. G. Taylor, “Neural modeling and functional brain imaging: an overview,” Neural Networks, vol. 13, no. 8-9, pp. 829–846, 2000.
View at: Publisher Site | Google Scholar
M. Ferrari and V. Quaresima, “A brief review on the history of human functional near-infrared spectroscopy (fNIRS) development and fields of application,” NeuroImage, vol. 63, no. 2, pp. 921–935, 2012.
View at: Publisher Site | Google Scholar
S. Debener, M. Ullsperger, M. Siegel, and A. K. Engel, “Single-trial EEG-fMRI reveals the dynamics of cognitive function,” Trends in Cognitive Sciences, vol. 10, no. 12, pp. 558–563, 2006.
View at: Publisher Site | Google Scholar
E. Martínez-Montes, P. A. Valdés-Sosa, F. Miwakeichi, R. I. Goldman, and M. S. Cohen, “Concurrent EEG/fMRI analysis by multiway partial least squares,” NeuroImage, vol. 22, no. 3, pp. 1023–1034, 2004.
View at: Publisher Site | Google Scholar
S. G. Horovitz, B. Rossion, P. Skudlarski, and J. C. Gore, “Parametric design and correlational analyses help integrating fMRI and electrophysiological data during face processing,” NeuroImage, vol. 22, no. 4, pp. 1587–1595, 2004.
View at: Publisher Site | Google Scholar
K. Mullinger, S. Debener, R. Coxon, and R. Bowtell, “Effects of simultaneous EEG recording on MRI data quality at 1.5, 3 and 7 tesla,” International Journal of Psychophysiology, vol. 67, no. 3, pp. 178–188, 2008.
View at: Publisher Site | Google Scholar
T. Yamamoto and T. Kato, “Paradoxical correlation between signal in functional magnetic resonance imaging and deoxygenated haemoglobin content in capillaries: a new theoretical explanation,” Physics in Medicine and Biology, vol. 47, no. 7, pp. 1121–1141, 2002.
View at: Publisher Site | Google Scholar
P. Rasmussen, E. A. Dawson, L. Nybo, J. J. van Lieshout, N. H. Secher, and A. Gjedde, “Capillary-oxygenation-level-dependent near-infrared spectrometry in frontal lobe of humans,” Journal of Cerebral Blood Flow & Metabolism, vol. 27, no. 5, pp. 1082–1093, 2007.
View at: Publisher Site | Google Scholar
G. Pfurtscheller, I. Daly, G. Bauernfeind, and G. R. Müller-Putz, “Coupling between intrinsic prefrontal Hbo2 and central EEG beta power oscillations in the resting brain,” PLoS ONE, vol. 7, no. 8, Article ID e43640, 2012.
View at: Publisher Site | Google Scholar
M. J. Herrmann, T. Huter, M. M. Plichta et al., “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 Mapping, vol. 29, no. 1, pp. 28–35, 2008.
View at: Publisher Site | Google Scholar
S. P. Koch, S. Koendgen, R. Bourayou, J. Steinbrink, and H. Obrig, “Individual alpha-frequency correlates with amplitude of visual evoked potential and hemodynamic response,” NeuroImage, vol. 41, no. 2, pp. 233–242, 2008.
View at: Publisher Site | Google Scholar
K. J. Worsley and K. J. Friston, “Analysis of fMRI time-series revisited—again,” NeuroImage, vol. 2, no. 3, pp. 173–181, 1995.
View at: Publisher Site | Google Scholar
T. Brookings, S. Ortigue, S. Grafton, and J. Carlson, “Using ICA and realistic BOLD models to obtain joint EEG/fMRI solutions to the problem of source localization,” NeuroImage, vol. 44, no. 2, pp. 411–420, 2009.
View at: Publisher Site | Google Scholar
G. H. Glover, “Deconvolution of impulse response in event-related BOLD fMRI,” NeuroImage, vol. 9, no. 4, pp. 416–429, 1999.
View at: Publisher Site | Google Scholar
K. J. Friston, A. P. Holmes, K. J. Worsley, J.-P. Poline, C. D. Frith, and R. S. J. Frackowiak, “Statistical parametric maps in functional imaging: a general linear approach,” Human Brain Mapping, vol. 2, no. 4, pp. 189–210, 1994.
View at: Publisher Site | Google Scholar
J. Cohen-Adad, S. Chapuisat, J. Doyon et al., “Activation detection in diffuse optical imaging by means of the general linear model,” Medical Image Analysis, vol. 11, no. 6, pp. 616–629, 2007.
View at: Publisher Site | Google Scholar
D. A. Handwerker, J. M. Ollinger, and M. D'Esposito, “Variation of BOLD hemodynamic responses across subjects and brain regions and their effects on statistical analyses,” NeuroImage, vol. 21, no. 4, pp. 1639–1651, 2004.
View at: Publisher Site | Google Scholar
G. M. Boynton, S. A. Engel, and D. J. Heeger, “Linear systems analysis of the fMRI signal,” NeuroImage, vol. 62, no. 2, pp. 975–984, 2012.
View at: Publisher Site | Google Scholar
T. J. Huppert, R. D. Hoge, S. G. Diamond, M. A. Franceschini, and D. A. Boas, “A temporal comparison of BOLD, ASL, and NIRS hemodynamic responses to motor stimuli in adult humans,” NeuroImage, vol. 29, no. 2, pp. 368–382, 2006.
View at: Publisher Site | Google Scholar
M. Takeuchi, E. Hori, K. Takamoto et al., “Brain cortical mapping by simultaneous recording of functional near infrared spectroscopy and electroencephalograms from the whole brain during right median nerve stimulation,” Brain Topography, vol. 22, no. 3, pp. 197–214, 2009.
View at: Publisher Site | Google Scholar
M. Tanosaki, Y. Hoshi, Y. Iguchi, Y. Oikawa, I. Oda, and M. Oda, “Variation of temporal characteristics in human cerebral hemodynamic responses to electric median nerve stimulation: a near-infrared spectroscopic study,” Neuroscience Letters, vol. 316, no. 2, pp. 75–78, 2001.
View at: Publisher Site | Google Scholar
P. Giacometti and S. G. Diamond, “Compliant head probe for positioning electroencephalography electrodes and near-infrared spectroscopy optodes,” Journal of Biomedical Optics, vol. 18, no. 2, Article ID 027005, 2013.
View at: Publisher Site | Google Scholar
Y. Hoshi, M. Shimada, C. Sato, and Y. Iguchi, “Reevaluation of near-infrared light propagation in the adult human head: implications for functional near-infrared spectroscopy,” Journal of Biomedical Optics, vol. 10, no. 6, Article ID 064032, 2005.
View at: Publisher Site | Google Scholar
A. J. Bell and T. J. Sejnowski, “A non-linear information maximisation algorithm that performs blind separation,” in Advances in Neural Information Processing Systems, vol. 7, pp. 467–474, MIT Press, 1995.
View at: Google Scholar
T.-P. Jung, S. Makeig, M. Westerfield, J. Townsend, E. Courchesne, and T. J. Sejnowski, “Removal of eye activity artifacts from visual event-related potentials in normal and clinical subjects,” Clinical Neurophysiology, vol. 111, no. 10, pp. 1745–1758, 2000.
View at: Publisher Site | Google Scholar
D. Mantini, M. G. Perrucci, S. Cugini, A. Ferretti, G. L. Romani, and C. del Gratta, “Complete artifact removal for EEG recorded during continuous fMRI using independent component analysis,” NeuroImage, vol. 34, no. 2, pp. 598–607, 2007.
View at: Publisher Site | Google Scholar
S. Hoffmann and M. Falkenstein, “The correction of eye blink artefacts in the EEG: a comparison of two prominent methods,” PLoS ONE, vol. 3, no. 8, Article ID e3004, 2008.
View at: Publisher Site | Google Scholar
P. Flandrin, G. Rilling, and P. Gonçalvés, “Empirical mode decomposition as a filter bank,” IEEE Signal Processing Letters, vol. 11, no. 2, pp. 112–114, 2004.
View at: Publisher Site | Google Scholar
S. Qian and D. Chen, “Joint time-frequency analysis,” IEEE Signal Processing Magazine, vol. 16, no. 2, pp. 52–67, 1999.
View at: Publisher Site | Google Scholar
P. S. Addison, “Wavelet transforms and the ECG: a review,” Physiological Measurement, vol. 26, no. 5, pp. R155–R199, 2005.
View at: Publisher Site | Google Scholar
H. Liang, S. L. Bressler, R. Desimone, and P. Fries, “Empirical mode decomposition: a method for analyzing neural data,” Neurocomputing, vol. 65-66, pp. 801–807, 2005.
View at: Publisher Site | Google Scholar
Z. Mehboob and H. Yin, “Information quantification of empirical mode decomposition and applications to field potentials,” International Journal of Neural Systems, vol. 21, no. 1, pp. 49–63, 2011.
View at: Publisher Site | Google Scholar
T. Talukdar, J. H. Moore, and S. G. Diamond, “Continuous correction of differential path length factor in near-infrared spectroscopy,” Journal of Biomedical Optics, vol. 18, no. 5, Article ID 056001, 2013.
View at: Publisher Site | Google Scholar
A. Villringer and B. Chance, “Non-invasive optical spectroscopy and imaging of human brain function,” Trends in Neurosciences, vol. 20, no. 10, pp. 435–442, 1997.
View at: Publisher Site | Google Scholar
J. P. Holman and W. J. Gajda, Experimental Methods for Engineers, vol. 2, McGraw-Hill, New York, NY, USA, 1994.
N. E. Huang and S. S. Shen, Hilbert-Huang Transform and Its Applications, vol. 5, World Scientific, 2005.
View at: MathSciNet
G. M. Boynton, S. A. Engel, G. H. Glover, and D. J. Heeger, “Linear systems analysis of functional magnetic resonance imaging in human v1,” The Journal of Neuroscience, vol. 16, no. 13, pp. 4207–4221, 1996.
View at: Google Scholar
M. A. Franceschini, S. Fantini, J. H. Thompson, J. P. Culver, and D. A. Boas, “Hemodynamic evoked response of the sensorimotor cortex measured noninvasively with near-infrared optical imaging,” Psychophysiology, vol. 40, no. 4, pp. 548–560, 2003.
View at: Publisher Site | Google Scholar
J. C. Lagarias, J. A. Reeds, M. H. Wright, and P. E. Wright, “Convergence properties of the Nelder—Mead simplex method in low dimensions,” SIAM Journal on Optimization, vol. 9, no. 1, pp. 112–147, 1998.
View at: Publisher Site | Google Scholar | MathSciNet
S. D. Muthukumaraswamy and B. W. Johnson, “Primary motor cortex activation during action observation revealed by wavelet analysis of the EEG,” Clinical Neurophysiology, vol. 115, no. 8, pp. 1760–1766, 2004.
View at: Publisher Site | Google Scholar
T. J. Huppert, S. G. Diamond, M. A. Franceschini, and D. A. Boas, “Homer: a review of time-series analysis methods for near-infrared spectroscopy of the brain,” Applied Optics, vol. 48, no. 10, pp. D280–D298, 2009.
View at: Publisher Site | Google Scholar
J. Han, M. Kamber, and J. Pei, Data Mining: Concepts and Techniques, Morgan Kaufmann, Boston, Mass, USA, 2006.
L. Ferreira and D. B. Hitchcock, “A comparison of hierarchical methods for clustering functional data,” Communications in Statistics. Simulation and Computation, vol. 38, no. 9, pp. 1925–1949, 2009.
View at: Publisher Site | Google Scholar | MathSciNet
R. V. Hogg and J. Ledolter, Applied Statistics for Engineers and Physical Scientists, vol. 59, Macmillan, New York, NY, USA, 1992.
G. Jasdzewski, G. Strangman, J. Wagner, K. K. Kwong, R. A. Poldrack, and D. A. Boas, “Differences in the hemodynamic response to event-related motor and visual paradigms as measured by near-infrared spectroscopy,” NeuroImage, vol. 20, no. 1, pp. 479–488, 2003.
View at: Publisher Site | Google Scholar
P. S. F. Bellgowan, Z. S. Saad, and P. A. Bandettini, “Understanding neural system dynamics through task modulation and measurement of functional MRI amplitude, latency, and width,” Proceedings of the National Academy of Sciences of the United States of America, vol. 100, no. 3, pp. 1415–1419, 2003.
View at: Publisher Site | Google Scholar
W. Richter, R. Somorjai, R. Summers et al., “Motor area activity during mental rotation studied by time-resolved single-trial fMRI,” Journal of Cognitive Neuroscience, vol. 12, no. 2, pp. 310–320, 2000.
View at: Publisher Site | Google Scholar
H. Girouard and C. Iadecola, “Neurovascular coupling in the normal brain and in hypertension, stroke, and Alzheimer disease,” Journal of Applied Physiology, vol. 100, no. 1, pp. 328–335, 2006.
View at: Publisher Site | Google Scholar
C. Iadecola, “Cerebral circulatory dysregulation in ischemia,” Cerebrovascular Diseases, vol. 1, pp. 319–332, 1998.
View at: Google Scholar

Copyright

Copyright © 2015 M. Tanveer Talukdar et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

1862

Downloads

1364

Citations