In this study a new strategy that combines Granger causality mapping (GCM) and independent component analysis (ICA) is proposed to reveal complex neural network dynamics underlying cognitive processes using functional near infrared spectroscopy (fNIRS) measurements. The GCM-ICA algorithm implements the following two procedures: (i) extraction of the region of interests (ROIs) of cortical activations by ICA, and (ii) estimation of the direct causal influences in local brain networks using Granger causality among voxels of ROIs. Our results show that the use of GCM in conjunction with ICA is able to effectively identify the directional brain network dynamics in time-frequency domain based on fNIRS recordings.
© 2013 Optical Society of America
Since the mid-1970, fNIRS has been employed for non-invasive investigation of brain activity by measuring the absorption coefficient of the near-infrared light between 650 nm and 950nm [1–9]. fNIRS offers unsurpassed high temporal resolution and provides quantitative information for both oxyhemoglobin (HbO2) and deoxyhemoglobin (HbR), which plays an important role in the in vivo study of cognitive processing in the human brain. Now advances in fNIRS are undergoing a transition from mapping sites of cortical activations towards identifying the brain networks that connect these sites together into dynamic systems in time-frequency domain. However, the algorithms for quantifying brain networks with fNIRS measurements are mainly confined to the temporal cross-correlation and the coherence spectrum analysis . Questions about revealing directional influences in anatomical and functional circuits are not addressed in these methods. Further, the time series data involved in the estimation of cross-correlation or coherence are typically limited to only two channels or bivariate analysis. It is known that these methods generally may not provide an entirely correct interpretation of inter-channel interactions, as they ignore influences from other channels [10,11].
Recently, a GCM method based on vector autoregressive models of fMRI/EEG time series was used to map effective connectivity in the human brain, promising further insights into the neural mechanisms of cognitive processing by revealing how one brain region might exert influence on another [10–12]. Directional information provided by Granger causality (G-C) offers the potential for defining the anatomical pathways that underlie neural interactions. Previous investigations have shown that GCM can resolve two important issues involved in the analysis of neural time series and network dynamics: (1) spectral analysis of fully multivariate (multichannel) time series and (2) spectral analysis of time series in a short time window . As such, GCM method has the capability for examining neurocognitive processing in brief time intervals to characterize its rapid temporal network dynamics and brain oscillations in frequency domain with multivariate analysis. However, the general GCM method is computationally expensive, as signals from thousands of voxels (or from all the data recording channels) have to be processed simultaneously in time-frequency domain to construct the brain networks and identify the brain dynamics at different time point and different frequency. In particular, GCM based on effective connectivity analysis method strongly depends on the accuracy of identified ROIs. If we miss the fNIRS neural sources or incorporate some unrelated physiology noise sources, the generated networks will not be reliable or effective.
Interestingly ICA has been recognized as a tool for evaluating the hidden spatiotemporal neural sources from EEG/fMRI measurements [12–14]. ICA also shows its potential for analyzing the independent sources to external stimuli in fNIRS recordings [5,6,9]. Basically ICA is a multivariate statistical method used to uncover brain activity sources from multiple data channels such that these sources are maximally independent. It is more interesting to look at time-frequency decompositions of component activations than of separate channel activity since independent component may directly index the activity of one fNIRS source of the brain, whereas channel activity acquires hemodynamic signals from different areas of the brain. As such, ICA is a very valuable tool to aid in the understanding of the brain cognitive processes in psychology or neurosciences .
In this study the GCM-ICA method is proposed for fNIRS data integration to identify brain network dynamics. The algorithm implements the following two procedures: (i) identification of the temporal time courses and the spatial activation maps (namely ROIs) from event-related brain activity by ICA and (ii) estimation of the direct causal influences in local brain networks using Granger causality among voxels/channels of ROIs. The developed fNIRS data integration method based on GCM-ICA will bring us a powerful tool to quantify the complex brain network dynamics in time-frequency domain with low computational cost and reliable accuracy, which provides the foundation for further investigations of neurological disorders and neural cognitive processes.
2.1 The Theory Model for fNIRS
According to Beer’s law , the wavelength-dependent tissue optical density changes can be written in terms of the changes of the chromophores including HbO2 and HbR at time t with wavelength λ,Eq. (1), the time series matrix for the changes of HbO2 and HbR is written as
ICA was first developed to tackle linear mixing problems similar to the ‘cocktail party problem’ in which many people are speaking at once and multiple microphones pick up different mixtures of the speakers' voices . The ICA algorithm employed in this study is infomax ICA algorithm, which addresses to separate mixed signals into maximally independent sources by maximization of information transfer between them [13–15].
Since fNIRS recordings are from K channels with M dimensional random time vector for each channel, the measurement matrix for is denoted by,12],
Step1: (e.g. random)
Step 3: If not converged, go back to step 2.in which t represents a given approximation step, is a general function that specifies the sizes of the steps for the unmixing matrix updates(usually an exponential function or a constant), I is the identity matrix, T is the transposition operator, and is the nonlinearity in the neural network.
In this study, we introduce the ICA method for fNIRS data processing. The significance of ICA is able to identify sites of cortical activations and characterize the hemodynamic task responses temporally and spatially for the whole brain, which establishes the foundation for further brain network dynamics analysis based on GCM method.
2.3 GCM Method
According to G-C, if the inclusion of past observations of X2 reduces the prediction error of X1 in a linear regression model of X1 and X2, as compared to a model which includes only previous observations of X1, X2 will cause X1. To illustrate G-C, suppose that the temporal dynamics of two time series X1(t) and X2(t) can be described by a bivariate autoregressive model :Eq. (8) (or Eq. (9)), we define X2 (or X1) causes X1 (or X2). Assuming that X1 and X2 are covariance stationarity, the magnitude of this interaction can be measured by the log ratio of the prediction error variances for the restricted (R) and unrestricted (U) models :Eq. (8) and E1U is derived from the full model in Eq. (8).
G-C can be easily extended to the multivariate(conditional) case in which the G-C of X2 on X1 is described in the context of multiple additional variables X3 . . . Xn. In this case, X2 causes X1 if knowing X2 reduces the variance in X1’s prediction error when all other variables X3 . . . Xn are also included in the regression model.
To illustrate, for a system of three variables (1,2,3), we represent the noise covariance matrix of the unrestricted model as
For multivariate autoregressive models, let be an l-dimensional random process (l is the voxels number from ROIs) defined in a segment of windowed time series data . Assuming stationarity of the process in each window, one can describe by a pth-order autoregressive process:
GCM method employs the Levinson-Wiggins-Robinson (LWR) algorithm to estimate the Am matrices and utilizes the Yule-Walker equations of the model to calculate the covariance matrix V of the noise vector E(t) . Then all the spectral quantities can be readily derived from the transfer functions using Fourier transformation:10]:
2.4 Behavior Tasks and fNIRS Recordings
The behavior test for fNIRS is implemented with a block design for a right finger tapping task. The experiment is performed using a 48-chnnel FOIRE-3000 (Shimadzu OMM) setup. As shown in Fig. 1, this system has 16 source, 16 detector, and 48 channels. In this system, two continuous wave lights at wavelengths 780nm and 856nm are emitted at each source fiber. The configurations of 48 channels within the cerebral cortex are provided in Fig. 2(a) for two-dimensional (2D) case and Fig. 2(b) for three-dimensional (3D) case, respectively. The raw fNIRS data sets are provided by Prof. Ye at KAIST (http://bisp.kaist.ac.kr/NIRS-SPM) and could be downloaded freely with permission.
In the case of block design for finger tapping tasks, the first onset time was right at 20s, then followed by a 20s period of activation alternated with a 40s period of rest. This was repeated 4 times for the subject. The last onset time was at 260s and the duration was 20s, then followed by a 20s period of rest. As such, the total recording time was 300s. The sampling rate of fNIRS recordings was 7.7Hz and DPF(r) = 4. During the task periods, subject was instructed to perform a finger flexion and extension action repeatedly.
Data segmentation, which is known as epoching in some data analysis systems, simply means chopping up the continuous fNIRS data into small time periods. The most common way to do this is to extract segments surrounding the event codes from the experiment (e.g., from −25s prior to the event code until 35s after the event code in this study). We adopt the method for fNIRS run epoching that has been extensively used in EEGLAB for EEG biopotential analysis (http://sccn.ucsd.edu/wiki/EEGLAB). To display the converted HbO2/HbR/HbT distributions, runs were imaged in (bottom-to-top) order of their occurrence during the experiment. For example, Fig. 3 shows their distributions for an arbitrary selected channel 36 that is close to the motor cortex .
ICA is an intrinsically multivariate approach, whose components provide a grouping of active brain regions that share the same response pattern. In this study the ROIs from fNIRS recordings are first extracted for right finger tapping tasks using ICA . Then the directed interactions among ROIs are explored using GCM analysis.
For the 4-run fNIRS recordings (5 onsets with period of 300s) of right finger tapping tasks, the 2D brain activation maps and their associated time courses were first identified by ICA from the raw HbO2 measurements (without filtering or pre-processing). To display the contribution of independent components (ICs) to the HbO2 recordings, we plot in Fig. 4(a) the components that contribute the most to HbO2 measurements in terms of the most HbO2 variance of all the 48 ICs. We also capture the components that contribute the most to the raw HbR and HbT measurements, which are provided in Fig. 4(b) and Fig. 4(c), respectively. It should be noted in Figs. 4(a)-4(c) that the black thick lines indicate the data envelope (i.e. minimum and maximum of all channels at every time point) and the colored ones show the HbO2/HbR/HbT components with the largest contributions to the measurements.
Then the most important ICs are sorted and grouped based on their contributions to generate the ROIs for each of the three chromophores. The identified ROIs by ICA consist of five nodes (voxels) for HbO2, six nodes for HbR, and five nodes for HbT, which are displayed on the fourth column of Figs. 5(a)-5(c), respectively. The nodes within the ROIs are considered as the primary cortex areas directly activated by the stimuli and will be employed for network analysis with GCM method. The VAR model order utilized in this study is 10, the window length are 80 time points and the frequency interval for fast Fourier transformation is between 1 and 100Hz. To investigate the brain networks underlying right finger tapping, the HbO2, HbR and HbT signals from all the nodes located within the ROIs are processed to calculate the time-frequency causal influences.
Directional influence measured as Granger causality using the time courses of raw HbO2 is provided in Fig. 6, in which Fig. 6(a) in the left column plots the time-frequency Granger causality influences among different nodes and Fig. 6 (b) in the right column shows the single spectra at 25 seconds (stimuli onset time). Likewise the computed time-frequency Granger causality influences for the HbR and HbT are displayed in Fig. 7 and Fig. 8, respectively.
As all node pairs show causal influences in both directions, a permutation procedure was applied to assign significance levels to the computed causalities. The significance thresholds corresponded to P < 0.05. Based on this, the Granger casual connectivity networks are reconstructed for arbitrarily selected time windows center at 25s (stimuli onset time) and 35s (stimulus processing time) for the results in Figs. 6-8, which are plotted in Fig. 9 for the three chromophores. The frequency selected for constructing the brain networks is 3 Hz, in which the brain activation patterns show significant correlations with the right finger tapping tasks, as displayed in Figs. 6-8. The edges of the brain networks are labeled with blue lines when they are estimated from HbO2/ HbR/ HbT signals, where the arrow represents the direction of Granger causality influences between two channels while the width of the blue lines represents the strength of Granger causality influences.
Mapping effective neuroconnectivity is an essential step toward unraveling the brain mechanisms of cognition and disorders. In this study, we have proposed a new approach combining GCM and ICA with fNIRS recordings to identify patterns of connectivity within a neural network implicated in right finger tapping processing.
The strength and direction of short-time causal influences were first examined in time-frequency domain during both the onset time and stimulus processing for all the node pairs within the ROIs. Significant causality influences in Figs. 6 and 8 are identified in the spectral vicinity of 3.0 Hz, with peaks generally occurring in the time period from 25s-45s (time points: 200-360), which are in good agreement with the stimulus duration of behavior tasks. However, this is not the case for HbR, in which strong causal influences are identified during the onset time around 25s while weak causality influences are revealed during the stimulus processing, as displayed in Fig. 7.
Importantly a number of features of causal influences are observed from the results plotted in Fig. 9. It is shown from the HbO2 networks in the first column of Fig. 9 that strong causal influences mainly occur in the left motor cortex (LMC) including Left primary motor cortex and left primary somatosensory cortex, which validate that stimuli in the right finger tapping tasks will yield the activations in LMC. In addition, we found in the first column of Fig. 9(a) that Granger causality influences are also observed from parietal cortex (PC) lobe region to LMC and from LMC to the supplementary motor area (SMA) during the onset time. These findings are consistent with the organization of motor cortex in its early stages. Interestingly during the stimulus processing, besides the strong Granger causality influences within LMC, increased causality influences among LMC, SMA and PC are observed from the second network of HbO2 at 35s, which demonstrates motor cortex information flows further to the motor association cortex such as PC. Further, we found during the stimulus processing SMA also exerts Granger causality influences on LMC. The generated effective networks for HbO2 correlate well with the structural networks and anatomical pathways in Fig. 10, in which the activation areas of motor stimuli are mainly located in the SMA and LMC for the right finger tapping tasks. In addition, brain activity is also observed in the motor association cortex such as PC during the stimulus processing, as displayed in Figs. 9 and 10. Interestingly we didn’t identify any neural activity from premotor cortex which generally involves the imaging of movement, but not the real motor stimuli such as finger tapping. The structural networks for finger tapping tasks could be found from the website (http://neuroscience.uth.tmc.edu/s3/chapter03.html).
It is seen from the HbR networks in the second column Fig. 9(a) that during the onset time strong causal influences mainly occur in LMC, but causal influences from LMC to SMA and PC are also identified. Moreover, we can see from the second column of Fig. 9(b) during the stimulus processing the influence within LMC is significantly reduced while influences from LMC to SMA are increased for HbR. In particular, the right primary motor cortex (RMC) also exerts influences on the LMC. These features are quite different from the findings in HbO2 networks.
Finally, if we compare the HbT networks in the third column of Fig. 9 with those from HbO2 and HbR, we found during both onset and stimulus processing, LMC exerts strong influences on SMC. And during the stimulus processing, strong Granger causality influence from LMC to PC is confirmed as well, as shown in the third column of Fig. 9(b).
We demonstrate in Figs. 6-8 how one can use a time-frequency analysis to potentially aid in the understanding of the neuromotor process underlying repetitive finger tapping. We found there is no neural oscillation activity over 4.5 Hz and synchronized delta frequency oscillations in the left primary motor cortex are involved in a simple finger tapping task for fNIRS recordings. In particular, the brain activations around 3Hz that correlate with stimulus processing are very significant compared those from 0 to 2Hz. The pattern of interaction we identified in frequency domain is in good agreement with previous fMRI studies, which validated the brain oscillations correlate well with the finger tapping frequency and are excellent features for the diagnosis of neurological disorders [12,17]. Importantly we observed from Figs. 6-8 the brain oscillation patterns are the same for the HbO2, HbR and HbT measurements though they show different activation patterns and network dynamics in temporal domain. It should be noted here that the ICs that reveal body movement and physiological noise such as ICs 16 and 44 in Fig. 11 are not involving the GCM analysis since only identified ROIs are employed for network analysis.
A key issue in effective connectivity mapping for fNIRS is how to define ROIs and how to represent the information in given ROIs in connectivity analysis. The classical GCM method is often employed to analyze of time series from all the voxels, which will result in high computational costs . However, the proposed approach first performs ICA analysis to extract ICs and identify ROIs. Since only signals from fNIRS channels/voxels within the ROIs, but not the whole channels/ voxels are processed, the computational cost could be reduced by over 90% compared with the general GCM method. As such combining ICA with GCM is very effective in detecting network connectivity while reducing the computational cost in fNIRS integration.
In conclusion, our new analysis method based on ICA and temporal-frequency Granger causality has shown to be useful in measuring complex effective connectivity from one brain region to the others with fNIRS recordings. The ICA method is crucial in enabling a whole brain network mapping with identified voxels by ICA at a reasonable computational cost. We thus suggest that the ICA-GCM technique is a potentially valuable tool that could be used for the investigation of causality influences of brain network dynamics in biophotonics fields.
This research was supported by MYRG Grant from University of Macau in Macau and Grant SRG2013-00035-FHS. In particular, we thank Prof. Mingzhou Ding, my previous colleague in Biomedical Engineering Department at University of Florida for providing me the Granger Causality method for my neuroscience projects. We also thank Prof. Jong Chui Ye at KAIST for providing us the fNIRS data sets for our initial methodology study.
References and links
2. B. Yodh and B. Chance, “Spectroscopy and imaging with diffusing light,” Phys. Today 48(3), 34–40 (1995). [CrossRef]
5. Z. Yuan, “Spatiotemporal and time-frequency analysis of functional near-infrared spectroscopy brain signals using independent component analysis,” J. Biomed. Opt. 18(10), 16011 (2013).
6. I. Schelkanova and V. Toronov, “Independent component analysis of broadband near-infrared spectroscopy data acquired on adult human head,” Biomed. Opt. Express 3(1), 64–74 (2012). [CrossRef] [PubMed]
7. R. C. Mesquita, M. A. Franceschini, and D. A. Boas, “Resting state functional connectivity of the whole head with near-infrared spectroscopy,” Biomed. Opt. Express 1(1), 324–336 (2010). [CrossRef] [PubMed]
8. 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,” Appl. Opt. 48(10), D280–D298 (2009). [CrossRef] [PubMed]
9. C. B. Akgül, A. Akin, and B. Sankur, “Extraction of cognitive activity-related waveforms from functional near-infrared spectroscopy signals,” Med. Biol. Eng. Comput. 44(11), 945–958 (2006). [CrossRef] [PubMed]
11. M. Ding, S. L. Bressler, W. Yang, and H. Liang, “Short-window spectral analysis of cortical event-related potentials by adaptive multivariate autoregressive modeling: data preprocessing, model validation, and variability assessment,” Biol. Cybern. 83(1), 35–45 (2000). [CrossRef] [PubMed]
12. O. Demirci, M. C. Stevens, N. C. Andreasen, A. Michael, J. Liu, T. White, G. D. Pearlson, V. P. Clark, and V. D. Calhoun, “Investigation of relationships between fMRI brain networks in the spectral domain using ICA and Granger causality reveals distinct differences between schizophrenia patients and healthy controls,” Neuroimage 46(2), 419–431 (2009). [CrossRef] [PubMed]
14. A. Delorme and S. Makeig, “EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis,” J. Neurosci. Methods 134(1), 9–21 (2004). [CrossRef] [PubMed]
15. Z. Yuan, J. Ye, “Fusion of fNIRS and fMRI data: Identifying when and where hemodynamic signal are changing in human brains,” Front. Hum. Neurosci., doi: 10.3389/ fnhum.2013.00676(2013).