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

Zhen Yuan

doi:10.1364/BOE.4.002629

1. 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–9]. fNIRS offers unsurpassed high temporal resolution and provides quantitative information for both oxyhemoglobin (HbO₂) 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–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–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.

2. Methods

2.1 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₂ and HbR at time t with wavelength λ,

[\begin{matrix} {Δ O D (r, t) |}_{λ_{1}} \\ {Δ O D (r, t) |}_{λ_{2}} \end{matrix}] = D P F (r) l (r) [\begin{matrix} ε_{1} (λ_{1}) & ε_{2} (λ_{1}) \\ ε_{1} (λ_{2}) & ε_{2} (λ_{2}) \end{matrix}] [\begin{matrix} Δ H b O_{2} (r, t) \\ Δ H b R (r, t) \end{matrix}]

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,

Δ O D

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,

Δ H b O_{2}

and

Δ H b R

(µ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₂ and HbR is written as

[\begin{matrix} Δ H b O_{2} (r, t) \\ Δ H b R (r, t) \end{matrix}] = [\begin{matrix} Q_{H b O_{2}} (r, t) \\ Q_{H b R} (r, t) \end{matrix}] / (D P F (r) l (r))

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:

[\begin{matrix} Δ H b O_{2} \\ Δ H b R \end{matrix}] = {(E^{T} R^{- 1} E)}^{- 1} E^{T} R^{- 1} [\begin{matrix} Δ O D |_{λ_{1}} \\ Δ O D |_{λ_{2}} \\ ... \\ Δ O D |_{λ_{n}} \end{matrix}] / (D P F \times l); E = [\begin{matrix} ε_{1} (λ_{1}) & ε_{2} (λ_{1}) \\ ... & ... \\ ε_{1} (λ_{n}) & ε_{2} (λ_{n}) \end{matrix}]

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

Δ H b T

(µM) is defined as the sum of

Δ H b O_{2}

and

Δ H b R

.

2.2 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 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 $Δ H b O_{2} (r, t)$ is denoted by,

Δ H b O_{2} = [\begin{matrix} Δ H b O_{2} (r_{1}, t_{1}) & Δ H b O_{2} (r_{1}, t_{2}) & ... & Δ H b O_{2} (r_{1}, t_{M}) \\ Δ H b O_{2} (r_{2}, t_{1}) & Δ H b O_{2} (r_{2}, t_{2}) & ... & Δ H b O_{2} (r_{2}, t_{M}) \\ ⋮ & ⋮ & ⋮ & ⋮ \\ Δ H b O_{2} (r_{K}, t_{1}) & Δ H b O_{2} (r_{K}, t_{2}) & ... & Δ H b O_{2} (r_{K}, t_{M}) \end{matrix}]

$Δ H b O_{2}$ 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,

Δ H b O_{2}

is further written [12],

Δ H b O_{2} = A S

in which A is the K-by-N mixing matrix, N is the number of unmixed sources and S are the N-by-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

W_{N \times K}

such that^{$X$ is given by $X = W (Δ H b O_{2} (r))$ 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: (e.g. random)
Step2:
$W (t + 1) = W (t) + η (t) (I - y (X) X^{T}) W (t)$ 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 X₂ reduces the prediction error of X₁ in a linear regression model of X₁ and X₂, as compared to a model which includes only previous observations of X₁, X₂ will cause X₁. To illustrate G-C, suppose that the temporal dynamics of two time series X₁(t) and X₂(t) can be described by a bivariate autoregressive model [16]:
$X_{1} (t) = \sum_{j = 1}^{p} A_{11, j} X_{1} (t - j) + \sum_{j = 1}^{p} A_{12, j} X_{2} (t - j) + E_{1} (t)$
$X_{2} (t) = \sum_{j = 1}^{p} A_{21, j} X_{1} (t - j) + \sum_{j = 1}^{p} A_{22, j} X_{2} (t - j) + E_{2} (t)$ in which p is the maximum number of lagged observations included in the model, A defines the coefficients of the model, and E₁ and E₂ are the variances for each time series. If the variance of E₁ (or E₂) is reduced by the inclusion of the X₂ (or X₁) terms in Eq. (8) (or Eq. (9)), we define X₂ (or X₁) causes X₁ (or X₂). Assuming that X₁ and X₂ 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]: $G_{2 \to 1} = \ln \frac{var (E_{1 R (12)})}{var (E_{1 U})}$ where E₁_R₍₁₂₎ is derived from the model omitting the A_12,_j (for all j) coefficients in Eq. (8) and E₁_U 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₂ on X₁ is described in the context of multiple additional variables X₃ . . . X_n. In this case, X₂ causes X₁ if knowing X₂ reduces the variance in X₁’s prediction error when all other variables X₃ . . . 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
$Σ = [\begin{matrix} var (E_{1 U}) & cov (E_{1 U}, E_{2 U}) & cov (E_{1 U}, E_{3 U}) \\ cov (E_{2 U}, E_{1 U}) & var (E_{2 U}) & cov (E_{2 U}, E_{3 U}) \\ cov (E_{3 U}, E_{1 U}) & cov (E_{3 U}, E_{2 U}) & var (E_{3 U}) \end{matrix}]$ where all E_i_U 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 $ρ = [\begin{matrix} var (E_{1 R}) & cov (E_{1 R}, E_{3 R}) \\ cov (E_{3 R}, E_{1 R}) & var (E_{3 R}) \end{matrix}]$ where all E_i_R 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
$G_{2 \to 1} |_{3} = \ln \frac{var (E_{1 R})}{var (E_{1 U})}$ For multivariate autoregressive models, let $X (t) = {[\begin{matrix} X_{1} (t) & X_{2} (t), & ...., & X_{l} (t) \end{matrix}]}^{T}$ be an l-dimensional random process (l is the voxels number from ROIs) defined in a segment of windowed time series data [11]. Assuming stationarity of the process $X (t)$ in each window, one can describe $X (t)$ by a pth-order autoregressive process:
$E (t) = \sum_{m = 0}^{p} A_{m} X (t - m)$ in which $A_{m},$ $m = [\begin{matrix} 1, & 2, & ..., & p \end{matrix}]$ are $l \times l$ coefficient matrix and $E (t)$ is a l-dimensional, zero mean and uncorrelated variance vector.
GCM method employs the Levinson-Wiggins-Robinson (LWR) algorithm to estimate the A_m matrices and utilizes the Yule-Walker equations of the model to calculate the covariance 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:
$H (f) = \sum_{m = 0}^{p} (A_{m} e^{- i m 2 π f})^{- 1}$ The spectral matrix of the time series data is given by $S (f) = H (f) V H^{*} (f)$ 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]:
$I_{j \to i} (f) = - \log (1 - \frac{(V_{j, j} - \frac{V_{i, j}^{2}}{V_{j, j}} {| H_{i, j} (f) |}^{2}}{S_{i, i} (f)})$ 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.
Fig. 1 Channel configurations of the Shimadzu’s FOIRE 3000 systems (Red circles: sources; blue circles: detectors).
Download Full
Size | PDF
Fig. 2 Channel configurations along the cerebral cortex for 2D and 3D analysis of fNIRS signals.
Download Full
Size | PDF
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₂/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].
Fig. 3 Run by run visualization in temporal domain for the raw HbO₂ (a), HbR (b) and HbT measurements (c) from an arbitrary selected channel 36. The bottom curve plots the run average for each chromophore. Runs were imaged in (bottom-to-top) order of their occurrence during the experiment.
Download Full
Size | PDF
3. 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₂ measurements (without filtering or pre-processing). To display the contribution of independent components (ICs) to the HbO₂ recordings, we plot in Fig. 4(a) the components that contribute the most to HbO₂ measurements in terms of the most HbO₂ 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₂/HbR/HbT components with the largest contributions to the measurements.
Fig. 4 The spatiotemporal analysis of independent components that contribute the most to the original HbO₂ measurements (a), HbR measurements(b), and HbT measurements (c). Brain activation maps in 2D are also provided and the cortical activations were mostly seen in the left primary motor cortex and supplementary motor area. (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 components of different chromophores).
Download Full
Size | PDF
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₂, 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₂, HbR and HbT signals from all the nodes located within the ROIs are processed to calculate the time-frequency causal influences.
Fig. 5 The most important independent components (ICs) calculated from HbO₂ measurements (a), HbR measurements (b), and HbT measurements (c). The combinations of these ICs for each specific chromophore will generate its ROIs for right finger tapping tasks and the ROIs are shown in the fourth column of Fig. 5. The cortical activations were mostly seen in the left primary motor cortex and supplementary motor area.
Download Full
Size | PDF
Directional influence measured as Granger causality using the time courses of raw HbO₂ 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.
Fig. 6 Directional influences measured as Granger causality among channels G1-G5 for HbO₂ measurements. Time-frequency causality influences are shown in the left column and the single spectra are given in the right column with the spectrum in the time window centered at 25s (time points: 200). Stimuli onset time: 25s (time point 200); stimulus processing time: 25s-45s (time points 200-360). Only significant influences are provided due to limited space.
Download Full
Size | PDF
Fig. 7 Directional influences measured as Granger causality among channels G1-G6 for HbR measurements. Time-frequency causality influences are shown in the left column and the single spectra are given in the right column with the spectrum in the time window centered at 25s. Stimuli onset time: 25s (time point 200); stimulus processing time: 25s-45 (time points 200-360). Only significant influences are provided due to limited space.
Download Full
Size | PDF
Fig. 8 Directional influences measured as Granger causality among channels G1-G5 for HbT measurements. Time-frequency causality influences are shown in the left column and the single spectra are given in the right column with the spectrum in the time window centered at 25s. Stimuli onset time: 25s (time point 200); stimulus processing time: 25s-45s(time points 200-360). Only significant influences are provided due to limited space.
Download Full
Size | PDF
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₂/ 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.
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₂ 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.
Download Full
Size | PDF
4. 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 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 HbO₂ 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₂ 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₂ 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).
Fig. 10 Motor cortex areas of human brain (open access website: http://neuroscience.uth.tmc.edu/s3/chapter03.html).
Download Full
Size | PDF
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₂ networks.
Finally, if we compare the HbT networks in the third column of Fig. 9 with those from HbO₂ 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 HbO₂, 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.
Fig. 11 The typical ICs that identify the body movement and boundary noise.
Download Full
Size | PDF
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.
Acknowledgments
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
1. F. F. Jöbsis, “Noninvasive, infrared monitoring of cerebral and myocardial oxygen sufficiency and circulatory parameters,” Science 198(4323), 1264–1267 (1977). [CrossRef] [PubMed]
2. B. Yodh and B. Chance, “Spectroscopy and imaging with diffusing light,” Phys. Today 48(3), 34–40 (1995). [CrossRef]
3. J. C. Ye, S. Tak, K. E. Jang, J. Jung, and J. Jang, “NIRS-SPM: statistical parametric mapping for near-infrared spectroscopy,” Neuroimage 44(2), 428–447 (2009). [CrossRef] [PubMed]
4. Z. Yuan, Q. Zhang, E. S. Sobel, and H. Jiang, “Image-guided optical spectroscopy in diagnosis of osteoarthritis: A clinical study,” Biomed. Opt. Express 1(1), 74–86 (2010). [CrossRef] [PubMed]
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]
10. J. Cui, L. Xu, S. L. Bressler, M. Z. Ding, and H. L. Liang, “BSMART: A MATLAB/C toolbox for analysis of multichannel neural time series,” Neural Netw. 21(8), 1094–1104 (2008). [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]
13. A. J. Bell and T. J. Sejnowski, “An information-maximization approach to blind separation and blind deconvolution,” Neural Comput. 7(6), 1129–1159 (1995). [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).
16. A. K. Seth, “A MATLAB toolbox for Granger causal connectivity analysis,” J. Neurosci. Methods 186(2), 262–273 (2010). [CrossRef] [PubMed]
17. W. Liu, L. Forrester, and J. Whitall, “A note on time-frequency analysis of finger tapping,” J. Mot. Behav. 38(1), 18–28 (2006). [CrossRef] [PubMed]}

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

Abstract

1. Introduction

2. Methods

2.1 The Theory Model for fNIRS

2.2 ICA

2.3 GCM Method

2.4 Behavior Tasks and fNIRS Recordings

3. Results

4. Discussion

Acknowledgments

References and links

Cited By

Figures (11)

Equations (17)

Biomedical Optics Express