Combining independent component analysis and Granger causality to investigate brain network dynamics with fNIRS measurements

: 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.


Introduction
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][2][3][4][5][6][7][8][9]. fNIRS offers unsurpassed high temporal resolution and provides quantitative information for both oxyhemoglobin (HbO 2 ) 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 [7]. 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][11][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 [10]. 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][13][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 [5].
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.

The Theory Model for fNIRS
According to Beer's law [1], the wavelength-dependent tissue optical density changes can be written in terms of the changes of the chromophores including HbO 2 and HbR at time t with wavelength λ, in which OD is the optical density as determined from the negative log ratio of the detected intensity of light with respect to the incident intensity of light using CW measurements, OD Δ is the optical density change (unitless quantity) at the position r, DPF(r) is the unitless differential path length factor, l(r) (mm) is the distance between the source and the detector, ( ) i ε λ is the extinction coefficient of the ith chromophore at wavelength λ of laser source, 2 HbO Δ and HbR Δ (µM) is the chromophore concentration change for oxy-and deoxyhemoglobin, respectively. After multiply the inverse matrix of the extinction coefficients for both sides of Eq. (1), the time series matrix for the changes of HbO 2 and HbR is written as HbO HbR Q r t HbO r t DPF r l r HbR r t Q r t in which Q (r,t) vectors are the product of the inversion matrix of the extinction coefficients and the optical density change vectors. Similar operational procedures could be extended to nth wavelength case based on regularization methods: where the matrix E is the extinction coefficient matrix and R is defined as the a priori estimate of the covariance of the measurement error. The change of total hemoglobin concentration HbT Δ (µM) is defined as the sum of 2 HbO Δ and HbR Δ .

ICA
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 [13]. The ICA algorithm employed in this study is ΔHbO is then decomposed by ICA, estimating the optimal inverse of the mixing matrix A, and a set of source time courses S. In terms of ICA estimation,

2
ΔHbO is further written [12], 2 ΔHbO = AS (5) in which A is the K-by-N mixing matrix, N is the number of unmixed sources and S are the Nby-M component time courses. Typically we utilize K ≥ N so that A is usually of full rank. The goal of ICA is to estimate an unmixing matrix N K W × such that X is given by 2 X = W(ΔHbO (r)) (6) in which X is a good approximation to the 'true' neural sources S and the inversion of weighted matrix W will extract the brain activity maps of the unmixed spatial sources. The infomax ICA algorithm uses the gradient ascent iteration algorithm to compute W by maximizing the entropy of the output of a single layer neural network. The resulting updated equation for the algorithm to calculate W is, Step1: Step2: Step 3: If not converged, go back to step 2.
in which t represents a given approximation step, ) (t η 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, 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.

GCM Method
According to G-C, if the inclusion of past observations of X 2 reduces the prediction error of X 1 in a linear regression model of X 1 and X 2 , as compared to a model which includes only previous observations of X 1 , X 2 will cause X 1 . To illustrate G-C, suppose that the temporal dynamics of two time series X 1 (t) and X 2 (t) can be described by a bivariate autoregressive model [16]: in which p is the maximum number of lagged observations included in the model, A defines the coefficients of the model, and E 1 and E 2 are the variances for each time series. If the variance of E 1 (or E 2 ) is reduced by the inclusion of the X 2 (or X 1 ) terms in Eq. (8) (or Eq. (9)), we define X 2 (or X 1 ) causes X 1 (or X 2 ). Assuming that X 1 and X 2 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 [16]: where E 1R (12) is derived from the model omitting the A 12,j (for all j) coefficients in Eq. (8) and E 1U 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 X 2 on X 1 is described in the context of multiple additional variables X 3 . . . X n . In this case, X 2 causes X 1 if knowing X 2 reduces the variance in X 1 's prediction error when all other variables X 3 . . . X n 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 where all E iU are estimated from the autoregressive model including all variables. For n variables there are n restricted models, with each restricted model omitting a different predictor variable. For example, the noise covariance matrix of the restricted model omitting variable 2, is where all E iR are estimated from the autoregressive model omitting variable 2. The G-C from variable 2 to variable 1, conditioned on variable 3, is given by For multivariate autoregressive models, let be an ldimensional random process (l is the voxels number from ROIs) defined in a segment of windowed time series data [11]. Assuming stationarity of the process ( ) t X in each window, one can describe ( ) t X by a pth-order autoregressive process: in which , matrix V of the noise vector E(t) [10]. Then all the spectral quantities can be readily derived from the transfer functions using Fourier transformation: The spectral matrix of the time series data is given by where * denotes matrix transposition and complex conjugation. Then the directional Granger casual influence from channel (voxel) j to channel i is given by for a specific frequency f [10]:

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 HbO 2 /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 [5].

Results
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 [5]. 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 HbO 2 measurements (without filtering or pre-processing). To display the contribution of independent components (ICs) to the HbO 2 recordings, we plot in Fig. 4(a) the components that contribute the most to HbO 2 measurements in terms of the most HbO 2 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 HbO 2 /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 HbO 2 , 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 HbO 2 , 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 HbO 2 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 HbO 2 / 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.

Discussion
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 timefrequency 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 HbO 2 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 HbO 2 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 HbO 2 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 HbO 2 networks.
Finally, if we compare the HbT networks in the third column of Fig. 9 with those from HbO 2 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). Fig. 9. Granger causality networks constructed for the time window centered at 25 s (a), and centered at 35s (b). 25s is the stimuli onset time and 35s is the stimulus processing time. The first column is for HbO 2 measurement, the second column is for HbR measurement while the third column is for HbT measurement. The arrow represents the direction of Granger causality influences while the width of the blue lines shows the strength of the influences. 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 HbO 2 , 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 [10]. 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.