2.1.

Theory Model for fNIRS

According to Beer’s law,¹ the wavelength-dependent tissue optical density changes can be written in terms of the concentration changes of the chromophores including ${\mathrm{HbO}}_2$ and HbR at time $t$ and wavelength $\lambda$,

2.2.

Independent Component Analysis

The ICA algorithm employed in this study is infomax ICA algorithm,¹⁸ which separates mixed signals into maximally independent sources by maximization of information transfer between them.¹⁹

2.3.

Experimental Systems and Data Acquisition

The fNIRS tests are implemented with a block design for a right finger tapping task. The experiment is performed using a 48-channel FOIRE-3000 (Shimadzu OMM) system, which has 16 sources, 16 detectors, and 48 channels. In this system, two continuous wave lights at wavelengths 780 and 856 nm are emitted at each source fiber. In the case of block design for right finger tapping tasks, the onset time for the first trigger was at 20 s and then followed by a 20-s period of activation alternated with a 40-s period of rest. This was repeated four times for the subject. The onset time for the last trigger was at 260 s and the duration was 20 s, followed by a 20-s period of rest. As such, the total recording time was 300 s. During the task period, the subject was instructed to perform a finger flexion and extension action repeatedly. Data segmentation, which is also known as epoching in signal processing, is utilized to chop up the continuous fNIRS data into small time periods. The general way to do this is to extract segments surrounding the event codes from the experiments, e.g., from $-25\mathrm{s}$ prior to the event onset until 35 s after the event code in this study. The original photon density datasets could be downloaded freely from http://bisp.kaist.ac.kr/NIRS-SPM/Sample_data.²⁰ The converted $\mathrm{Hb}{\mathrm{O}}_2$, HbR, and HbT change measurements ($\mathrm{DPF}=4$, $\text{sampling rate}=7.7\mathrm{Hz}$) are filtered, segmented, and displayed in Figs. 1(a)–1(c), respectively. The configurations of 48 channels located on the scalp are also provided in Fig. 1(d). It is seen from Fig. 1 that significantly higher increases in $\mathrm{\Delta }\mathrm{HbT}$ and $\mathrm{\Delta }{\mathrm{HbO}}_2$ in the motor cortex were observed during the stimulus processing for the finger tapping tasks in comparison with those from rest states. However, this is not the case for $\mathrm{\Delta }\mathrm{HbR}$, in which we found a significant decrease during the stimulus processing when compared with the measurement from baselines.

Fig. 1

The data reviews for filtered $\mathrm{Hb}{\mathrm{O}}_2$ measurements (a), HbR measurements (b) and HbT measurements (c). Channel configurations in two-dimension (2-D) (left column) and three-dimension (3-D) (right column) (d). Here, four-run fNIRS recordings (separated by dashed lines) are plotted at 48 channel sites (channel names on the left) and the onset time for right finger tapping is marked as “target.”

3.1.

Spatiotemporal Analysis of fNIRS Signals Using ICA

For all the four-run fNIRS recordings (five onsets with period of 300 s) of right finger tapping tasks, the 2-D brain activity maps and associated IC time courses were first identified by ICA from the filtered $\mathrm{Hb}{\mathrm{O}}_2$, HbR and HbT change measurements. To quantify the contribution of ICs to $\mathrm{Hb}{\mathrm{O}}_2$, we plot in Fig. 2(a) the components that contribute the most to it in terms of the most $\mathrm{Hb}{\mathrm{O}}_2$ variance of all the 48 ICs, where we observed ICs 1–4, 6, 7, and 14 are the most significant ones. However, we found that only ICs 1, 4, 6, and 14 correlate well with right finger tapping tasks since ICs 2, 3, and 7 basically identify noise or unrelated neural events. The brain activity maps of ICs 1, 4, 6, and 14 are the identified ROIs for $\mathrm{Hb}{\mathrm{O}}_2$. We also capture the components that contribute the most to the filtered HbR and HbT measurements, which are provided in Figs. 2(b) and 2(c), respectively.

Fig. 2

The spatiotemporal analysis of independent components (ICs) that contribute the most to the $\mathrm{Hb}{\mathrm{O}}_2$ measurements (a), HbR measurements (b), and HbT measurements (c). Brain activity maps in 2-D are also provided and the cortical activity was mostly seen in the left primary motor cortex and supplementary motor area (SMA) (see detailed images in Fig. 3). The black thick lines indicate the data envelope (i.e., minimum and maximum of all channels at every time point) and the colored lines show the components of different chromophores.

In addition, the ROIs identified by ICA from different chromophores are provided in Fig. 3, where we found that the ROIs for right finger tapping tasks are mainly located in the left primary motor cortex and supplementary motor area. However, we did observe identified brain activities located in the parietal cortex for $\mathrm{Hb}{\mathrm{O}}_2$ or in the right motor cortex for HbR. These findings show good agreement with neural physiologies for a right finger tapping task.^3,9 We also reconstructed the brain images for a different subject using fMRI data with the right finger tapping tasks (fMRI datasets can be downloaded from http://bisp.kaist.ac.kr/NIRS-SPM).²¹ We found the ROIs of fNIRS correlate very well with fMRI findings while the recordings are acquired from different subjects, which validate that ICA is able to effectively capture the spatial maps of brain activity using fNIRS recordings. Further work that combines the IC time courses and NIRS-SPM⁹ is being done to achieve high resolution and 3-D ICs imaging.

Fig. 3

The most significant ICs calculated from $\mathrm{Hb}{\mathrm{O}}_2$ measurements (a), HbR measurements (b), and HbT measurements (c). The composite ICs for each specific chromophore will generate its ROIs for right finger tapping tasks. (d) The fMRI image for one subject who performed right finger tapping tasks. The cortical activity was mostly seen in the left primary motor cortex and SMA and the fNIRS findings correlated well with the fMRI features. LM: left primary motor cortex; LM1: left primary motor cortex 1; LM2: left primary motor cortex 2; LM3: left primary motor cortex 3; SMA: supplementary motor area; RM: right primary motor cortex; P: parietal cortex; and V: visual cortex.

To display the hemodynamic response distributions processed by ICA, runs were imaged in (bottom-to-top) order of their occurrence during the experiment. Figure 4 shows their distribution in an arbitrary selected channel 41 after we remove the unrelated ICs including artifacts or physiological noise. The bottom curve of each figure in Fig. 4 shows the run-averaged time series of fNIRS recordings for each chromophore. Compared with the datasets processed without ICA shown on the right column of Fig. 4, the new datasets on the left column show correct run duration (0 to 20 s) and accurate onset time (0 s), which validates that the ICA has the capabilities for extracting the correct components during stimulus processing in temporal domain.

Fig. 4

Run by run visualization in temporal domain for the $\mathrm{Hb}{\mathrm{O}}_2$ measurements (a), HbR measurements (b) and HbT measurements (c) processed with ICA (left column of the images) and without ICA (right column of the images). The images on the left column show the correct onset time (0 s) and stimulus duration (0 to 20 s). Runs were imaged in (bottom-to-top) order of their occurrence during the experiment for an arbitrary selected channel 41. The bottom curve of each figure shows the run-averaged time series of fNIRS recordings for each chromophore.

3.2.

Time-Frequency Analysis of fNIRS Signals Using ICA

It is more interesting to look at time-frequency decompositions of component activations than of separate channel activity since ICs may directly index the activity of one brain fNIRS source, whereas channel activities acquire hemodynamic signals from different areas of the brain. We demonstrate how one can use a time-frequency domain ICA to potentially aid in the understanding of the neuromotor process underlying repetitive finger tapping. In Figs. 5(a)–5(c), we plot the ICs that account for the largest portions of brain oscillations with frequencies at 3.5 and 1.2 Hz for $\mathrm{Hb}{\mathrm{O}}_2$, HbR, and HbT, respectively. The ICs in frequency domain are sorted and selected in terms of their percentage contributions of total data power. Our studies have shown that synchronized delta frequency (0 to 3 Hz) oscillations are involved in a simple finger tapping task with fNIRS measurements. In particular, the brain oscillatory activities around $1.7/2.6/3.5\mathrm{Hz}$ are very significant compared those from other frequency bands. Interestingly, it was observed from the right column of Fig. 5 that the ROIs for right finger tapping tasks around 3.5 Hz are found mainly in the left primary motor cortex. However, this is not the case for the brain activation patterns around 1.2 Hz, in which ROIs are not related with the left motor cortex, as displayed on the left column of Fig. 5. It is observed from Fig. 5 that the most significant brain oscillatory activity was mainly found in the left primary motor cortex, which is well correlated with the right finger tapping tasks.

Fig. 5

The ICs accounting for the largest portions of powers of brain oscillatory activities with frequencies at 1.2 Hz (images on the left column) and 3.5 Hz (images on the right column) for $\mathrm{Hb}{\mathrm{O}}_2$ measurements (a), HbR measurements (b) and HbT measurements (c). The brain oscillations at 3.5 Hz correlated very well with the right finger tapping tasks for all of the three chromophores. LM: left primary motor cortex and SMA: supplementary motor area.

ERSP is able to capture the mean changes in spectral power at each time during the run and at each frequency. In terms of the ERSP distributions shown in Fig. 6, representative IC4 (the brain activity from the fourth IC is well correlated with right finger tapping tasks) from $\mathrm{Hb}{\mathrm{O}}_2$ shows increased ERSP near delta frequency band (1 to 3 Hz) during the stimulus processing (5 to 30 s), whereas representative IC11 from HbR shows a power decrease in the time range ($-5$ to 10 s), then followed by an increase during the stimulus processing (10 to 30 s). According to the distributions of ITCs plotted in Fig. 6 the brain activities from both $\mathrm{Hb}{\mathrm{O}}_2$ and HbR components appear to become partially synchronized around 3.5 Hz during the stimuli period from 0 to 15 ms after the onset of finger tapping tasks. However, brain activities for IC4 from $\mathrm{Hb}{\mathrm{O}}_2$ become partially synchronized around 1.7 Hz during the stimulus period (15 to 30 s) while IC11 from HbR becomes partially synchronized around 1.7 Hz before the onset of trigger ($-15$ to 0 s) and around $1.7/2.6\mathrm{Hz}$ during the stimulus processing.

Fig. 6

Event-related spectral perturbation (ERSP) plots the mean changes in spectral power at each time during the run and at each frequency. The IC4 (the fourth IC) of $\mathrm{Hb}{\mathrm{O}}_2$ shows increased ERSP during the stimulus processing with frequencies at $0.7/1.7/2.6/3.5\mathrm{Hz}$, whereas IC11 of HbR first shows the decreased ERSP centered at the onset time and then later increased ERSP till the end of the duration with frequencies at $0.7/1.7/2.6/3.5\mathrm{Hz}$. The brain activities for IC4 and IC11 become partially synchronized around 3.5 Hz during the stimulus period of 0 to 15 s. However, the brain activities for IC4 become partially synchronized around 1.7 Hz during the stimulus period of 15 to 30 s, whereas IC11 from HbR becomes partially synchronized around 1.7 Hz before the onset of trigger ($-15$ to 0 s) and around $1.7/2.6\mathrm{Hz}$ during the stimulus processing. The left curve in the ERSP panel shows the baseline mean power spectrum, whereas the lower curve of the ERSP panel shows the ERSP envelope [low (in blue) and high (in green) mean decibel values, relative to baseline, at each time in the run]. The left curve in intertrial coherence (ITC) panel shows the distributions of mean ITC at different frequencies, whereas the bottom curve in ITC panel shows the run-averaged time course of IC4 or IC11 in micormolar. (a) The computed ERSP and ICTC for IC4 from $\mathrm{Hb}{\mathrm{O}}_2$ measurements and (b) the computed ERSP and ITC for IC11 from HbR measurements.

3.3.

Comparisons of Different ICA Methods for fNIRS Signal Analysis

The MATLAB tool packages have over 20 available ICA algorithms.²² We calculated and compared the recovered ICs using four popular ICA methods including infomax, jade, sobi, and acsobiro. The calculated results using infomax based on HbT measurements are shown in Fig. 7(a), whereas Figs. 7(b)–7(d) show the most important ICs computed from methods 2 to 4. By estimating the spatiotemporal profiles of ICs, we found the recovered ICs using infomax are in good agreement with those from methods 2 to 4. The brain activity areas identified by methods 1 to 4 were mainly found in the left motor cortex, which correlate well with the neurophysiology mechanism for a right finger tapping task. Consequently, there are no significant differences for the reconstructed ICs from the four methods though minor calculation errors among them do exist. It seems that infomax, jade, and acsobiro capture similar distribution of HbT components compared with that from the sobi method. If we compare the recovered ICs 2 and 5 from jade with similar ICs 1 and 5 from infomax, ICs 3 and 1 from sobi, and ICs 1 and 5 from acsobiro, we would find jade identifies more features than the other methods. In addition, we also observed that infomax is able to extract the spatial activity maps with high contrast though with strong background noise. The main advantages of infomax algorithms based on the minimization of mutual information are its ability to adapt to variations in the environment and the fact that it is robust if the right type of distribution is provided (super- or sub-Gaussian).²⁰

Fig. 7

The spatiotemporal analysis of ICs that contribute the most to the HbT measurements using infomax (a), jade (b), sobi (c), and acsobiro (d) ICA methods. The black thick lines indicate the data envelope (i.e., minimum and maximum of all channels at every time point) and the colored lines show the components of HbT. No significant differences for the calculated ICs among the four methods were observed.

Both principal component analysis (PCA) and ICA could be used to extract independent neural sources or reduce the measurements into a smaller set of components. There are significant differences between the two blind source separation methods: (1) PCA uses the first and second moments of the measurements, hence relying heavily on Gaussian features while ICA exploits inherently non-Gaussian features of the data and utilizes higher order moments; (2) PCA maximizes the variance of the projected data along orthogonal direction while ICA finds the vectors onto which the projections are independent; and (3) in PCA, the first principle component accounts for as much of the variability in the data as possible, and each successive orthogonal component accounts for as much of the residual variability as possible while with ICA, we must first choose the number of sources to compute. Both ICA and PCA could be applied to the same problem and the results would be quite different. For example, in “cocktail party effect,”^19,20 ICA is able to distinguish the voice of each independent speaker from the linear combination of the voices while PCA fails to do that. Importantly, recent work¹³ on fNIRS shows that ICA has significant advantages on source separations in comparison with PCA.

It is noted that five conditions should be satisfied to perform ICA and PCA: the source signals must be statistically independent; the number of source signals must be equal to the number of mixed observed signals and mixtures must be linearly independent from each other; the model must be noise free; data must be centered; and the source signals must not have Gaussian probability density functions. Further investigations on ICA should include ICA for nonlinear mixing process, ICA for source signals that are noisy, ICA for a number of source signals greater than the number of observables and blind source separation techniques based on temporal dependencies.²⁰

In conclusion, we introduce ICA for fNIRS brain signal processing. The significance of this work is that ICA is able to identify sites of cortical activations and characterize the hemodynamic task responses in time-frequency domain for the whole brain with the event-related components. With further application to more complex tasks, ICA will likely reveal brain dynamics not identified with conventional fNIRS analysis methods.