2.1.

Basis of the Modified Beer-Lambert Law of Functional Near-Infrared Spectroscopy

From the modified Beer-Lambert law, the following equation is derived:

For multiple measurements at three wavelengths (780, 805, and $830\mathrm{nm}$ ), the following equation is derived from Eq. 2:

Then, with the extinction coefficients reported by Matcher, ² the following equation is derived using the least-squares method:

2.2.

Functional Near-Infrared Spectroscopic Imaging System

We used a 36-channel fNIRS imaging system with 12 source and 12 detector optodes (OMM-2001, Shimadzu Company, Kyoto, Japan).^{18, 19} The illuminator bundle quartz fiber (diameter $1\mathrm{mm}$ ) of the fNIRS imaging system was connected to three laser diodes (780, 805, and $830\mathrm{nm}$ ). After $10\mathrm{ms}$ of dark signal acquisition, each laser diode sequentially emitted a 5-ms pulse of light, and the light emission then traveled to a source optode. The sequence of light emission from the 12 source optodes (T1 to T12) was programmable.

2.3.

Flexible Adjustable Surface Holder

Holder caps made from thermoplastic resin are typically used to hold the optodes of an fNIRS imaging system, because the distance between the source and detector optodes must be constant. However, not all of our subjects could properly fit a cap onto their heads because of size differences among the individuals. This caused deterioration in the signal-to-noise ratio (SNR) and lowered the data reliability. We therefore developed a flexible adjustable surface holder (FLASH) to solve these problems [Fig. 1a ].

Fig. 1

Arrangement of fNIRS optodes and locomotor tasks. (a) Arrangement of source and detector optodes with the flexible adjustable surface holder (FLASH). (b) Locomotor tasks on the treadmill. Changes in hemodynamic response in the front-parietal area of healthy subjects were simultaneously recorded using an fNIRS imaging system (OMM-2001) and a laser Doppler tissue blood flow meter (FLO-C1). (c) Configuration of source and detector optodes in this study. (d) 24 optodes (12 source and 12 detector optodes) that enabled 36-channel recording were applied to the scalp of the front-parietal area.

The basic structure of the FLASH consists of sides and nodes making up triangles or quadrangles. Optode sockets are located on the nodes. The material of the sides is flexible, but not stretchable. The sum of the interior angles can be changed by allowing deflection of the sides and their rotation at the socket positions. Therefore, the FLASH shape can be defined by setting the interior angles at the sockets based on the Gauss-Bonnet theorem. We could thus change the FLASH shape and size to fit the various head shapes of our subjects without changing the distance between the source and detector.

2.4.

Theory of the Laser Doppler Tissue Blood Flow Meter

A laser Doppler tissue blood flow meter (FLO-C1 Omegawave Incorporated, Tokyo, Japan) measures tissue blood flow continuously and noninvasively at a detection depth of approximately $1\mathrm{mm}$ from the tissue surface. The tissue is irradiated by laser light $\left(780\mathrm{nm}\right)$ delivered through an optical fiber, and the light is scattered many times in the tissue. Some of the photons are scattered by erythrocytes in the tissue, and the light frequency shifts depending on the velocity of the erythrocyte movement. In addition, the amplitude of the scattered light is proportional to the volume of erythrocytes in the tissue. Blood flow signals are obtained from the frequency shift and its amplitude.

2.5.

Measurements

We chose to use a locomotor task for this study, because skin blood flow is changed by thermal skin stimulation. To answer our first question, we designed a 9-km/h locomotor task on a treadmill [Fig. 1b] that was expected to evoke large changes in the skin blood flow signal. The skin blood flows in the foreheads of healthy human subjects were measured using the laser Doppler tissue blood flow meter. A block design consisting of prerest $\left(20\mathrm{s}\right)$ , task $\left(20\mathrm{s}\right)$ , and postrest $\left(20\mathrm{s}\right)$ with three repetitions was used as the protocol. Six right-handed, healthy males who were 37 to 52 years of age participated in the study. The sensor of the laser tissue blood flow meter was located at the center of the forehead, $20\mathrm{mm}$ above the eyebrow line. We set the sampling rate of the laser Doppler tissue blood flow meter to $100\mathrm{ms}$ .

The 9-km/h speed chosen for the locomotor task was an experimental maneuver to get reproducible results for our first question. Though a paper using a 9-km/h locomotor task has been published,¹⁹ little knowledge about the brain signal was obtained compared with the use of a 3-km/h task. In addition, the task was excessive and susceptible to motion artifacts. Therefore, to answer our second question, during a 3-km/h locomotor task on a treadmill, changes in the hemodynamic response in the front-parietal area of two right-handed, healthy males were simultaneously recorded using the fNIRS imaging system and the laser Doppler tissue blood flow meter. ICA was applied to the fNIRS signals. A block design consisting of prerest $\left(30\mathrm{s}\right)$ , task $\left(90\mathrm{s}\right)$ , and postrest $\left(30\mathrm{s}\right)$ with three repetitions was used as the protocol. 24 optodes (12 source and 12 detector optodes), enabling 36-channel recording, were applied to each subject’s scalp [Figs. 1c and 1d]. The source optode (T8) next to the posterior one in the center row was located in the Cz portion based on the international 10 to 20 system for electroencephalogram electrode placement.⁵³ The interoptode distance was set to $30\mathrm{mm}$ . The sensor of the laser tissue blood flow meter was located at the center of the forehead, $20\mathrm{mm}$ above the eyebrow line. To reveal the influence of white noise, it was input into channels 30 and 36, because they were not located within an area of interest during the locomotor task. We acquired data by sequentially emitting laser light from T1 to T12 in numerical order. Thus, the fNIRS time resolution was $190\mathrm{ms}$ (see Sec. 2.2). In addition, we set the sampling rate of the laser Doppler tissue blood flow meter to $200\mathrm{ms}$ , because the fNIRS time resolution was $190\mathrm{ms}$ .

The study protocol was approved by the local ethics committee. Before the experiments, written informed consent was obtained from the subjects after they had received a complete description of the study.

2.6.

Independent Component Analysis

ICA assumes that the observed data vector $\mathbf{x}\left(t\right)=\left[{\mathbf{x}}_1\right(t),{\mathbf{x}}_2(t),\dots ,{\mathbf{x}}_n(t){]}^T:(t=1,2,\dots ,p)$ is a linear combination of unknown and statistically independent sources $\mathbf{s}\left(t\right)=\left[{\mathbf{s}}_1\right(t),{\mathbf{s}}_2(t),\dots ,{\mathbf{s}}_n(t){]}^T$ . The mixing matrix describing the linear combination of $\mathbf{s}\left(t\right)$ is given by $\mathbf{A}$ , such that

The algorithms must find a separating or demixing matrix $\mathbf{W}$ , such that

The ICA algorithm proposed by Molgedey and Schuster (MS-ICA)⁴⁰ decomposes $n$ -channel data into $n$ independent components, each of which corresponds to a recovered putative source that contributes to fNIRS signals. The main advantage of MS-ICA is its noniterative nature, and hence its high speed, but it is not directly guaranteed to minimize the likelihood. MS-ICA exploits the time coherence of the source signals to decompose the mixture of sources. It finds $\mathbf{W}$ by minimizing the sum-squared cross-correlations between one component at time $t$ and another component at time $t+\tau$ , across a set of time delays.⁴³ For various time delays, MS-ICA has been shown to be robust when the SNR is relatively poor.⁵⁴

We applied MS-ICA to the 36 channels of fNIRS signals. The ICA algorithm used here was implemented by ourselves. In the present study, the following time delays were tested: $\tau =0.19$ , 0.38, 0.57, 0.76, 0.95, 1.04, 1.33, 1.52, 1.71, 1.90, 2.19, 2.38, 2.57, 2.76, 2.95, 3.04, 3.33, 3.52, 3.71, and 3.80. They were chosen to cover a reasonably wide interval without extending beyond the support of the autocorrelation function.

3.1.

Skin Blood Flow

When we measured the skin blood flow in the forehead of six right-handed, healthy males during a 9-km/h locomotor task with a laser Doppler tissue blood flow meter (Fig. 2 ), we clearly observed changes in the skin blood flow signal depending on the task being performed for all subjects ( $\mathrm{p}<0.000001$ for all subjects). The statistical evaluation was based on a least-squares estimation using the general linear model. The design matrix was generated with a box-car function convolved with a Gaussian kernel of dispersion with a full width at half maximum of $4\mathrm{s}$ .

Fig. 2

One representative example of changes in the skin blood flow during a 9-km/h locomotor task on a treadmill. The skin blood flow on the forehead of healthy human subjects was monitored with a laser Doppler tissue blood flow meter. The skin blood signal changes depending on the task performed were observed for all subjects ( $\mathrm{p}<0.000001$ for all subjects).

When we measured the skin blood flow in the forehead of two right-handed, healthy males during a 3-km/h locomotor task, we observed changes in skin blood flow signal depending on the task being performed. The skin blood flow signal of one subject showed large increases with almost the same amplitude during all three task periods. However, the skin blood flow signal of another subject showed small increases during the first and second task periods, but a large increase during the third task period (Fig. 3 ). These results suggest that the skin blood flow was not always directly related to the locomotor task, but was indirectly related to the task via the autonomic nerve system. This observation supports the idea that changes in skin blood flow are caused by changes in skin temperature and sweating.^{32, 33, 34}

Fig. 3

One example of changes in skin blood flow during a 3-km/h locomotor task. The blood flow on the forehead of a healthy subject was monitored with a laser Doppler tissue blood flow meter. In this example, the skin blood flow signal showed small increases during the first and second task periods, but a large increase during the third task period.

3.2.

Independent Component Signals of the Functional Near-Infrared Spectroscopic Signals

During brain activation, an increase in the regional cerebral blood flow accompanies an increase in oxygenated hemoglobin. ^{4, 5, 6, 7, 8} Therefore, we analyzed the time course of $\mathrm{\Delta }\mathrm{oxy}$ calculated from fNIRS data measured simultaneously with the skin blood flow shown in Fig. 3. Figure 4a shows a time course of $\mathrm{\Delta }\mathrm{oxy}$ of a representative channel [channel 23 in Figs. 1c and 1d] corresponding to a foot area within the primary motor cortex. Figure 4b shows the power spectrum in the frequency range from 0 to $0.2\mathrm{Hz}$ in Fig. 4a calculated by Fourier transform. Figure 4c shows the power spectrum in the frequency range from 0.2 to $2\mathrm{Hz}$ in Fig. 4a calculated by Fourier transform. A large peak was found at around $0.006\mathrm{Hz}$ , which corresponded to the frequency calculated from our locomotor task period (i.e., $1∕150\mathrm{s}$ ), and would express the changes in skin blood flow as well as brain activation, depending on the task being performed. In addition, some peaks were found at around $0.08\mathrm{Hz}$ [low frequency oscillations, (LFO)] and $0.04\mathrm{Hz}$ [very low frequency oscillation, (VLFO)].²⁹ The cardiac pulsation frequency in our locomotor task was found to be in the frequency range from 1.5 to $1.8\mathrm{Hz}$ , but the respiration frequency was not found to be in the frequency range from 0.2 to $0.5\mathrm{Hz}$ . It is reasonable to hypothesize that changes in the skin blood flow and the regional cerebral blood signal are statistically independent, because these innervations are completely different. Hence, if the fNIRS signals included both changes in skin blood flow and cerebral blood signal, the signal changes from the two blood signals should have been separated by ICA. Even if changes in skin blood flow were dependent on the cerebral blood signal in our experimental coordinate system, ICA would have caused their interference to vanish in the new coordinate system generated ICA, because we hypothesize that they are statistically independent. The validity of the hypothesis will be checked by the consistency of the results. To remove various noises in fNIRS signals, it might be effective to use some processing such as filtering before ICA. However, the use of data processing other than ICA would not clarify the logical consequence of ICA. The data processing might distort the fNIRS signals. Therefore, only the MS-ICA algorithm, which is robust against noise, was used without filtering, smoothing, or baseline correction before the ICA. We processed the time courses of $\Delta \left[\mathrm{oxy}\text{-}\mathrm{Hb}\right]$ of the 36 channels using MS-ICA. The relevant color maps of a mixing matrix and the signal corresponding to independent components separated by the ICA are shown in Figs. 5 and 6 , respectively. Each color map in Fig. 5 represents the activated spatial distribution of the corresponding component, and the distribution corresponds to the arrangement of fNIRS channels shown in Figs. 1c and 1d. Note that they are spatial patterns without temporal information. The activated spatial distributions of components 1 and 3 in Fig. 5 were localized in the areas corresponding to channels 30 and 36, into which white noise was input. The results of spectral analysis for the corresponding components in Fig. 6 did not reveal any peaks. Therefore, we speculate that the components of the white noise were separated by ICA. The activated spatial distributions of components 2, 4, 6, 7, 9, and 11 in Fig. 5 show some characteristic patterns. The results of spectral analysis for the corresponding components in Fig. 6 had one or two peaks around $0.01\mathrm{Hz}$ (data not shown). Some of these components seem to be task dependent. They may be interpreted as some independent components by the difference in brain function. However, in the present study, we focused on determining the component related to changes in skin blood flow to remove it as an artifact instead of investigating all the components.

Fig. 4

Original time course of $\mathrm{\Delta }\mathrm{oxy}$ [defined by Eq. 4] and its power spectrum. (a) Time course of $\mathrm{\Delta }\mathrm{oxy}$ of a representative channel [channel 23 shown in Figs. 1c and 1d] corresponding to a foot area within the primary motor cortex. (b) The power spectrum (0 to $0.2\mathrm{Hz}$ ) of Fig. 4a calculated by Fourier transform. A large peak was found at around $0.006\mathrm{Hz}$ , which is the frequency calculated from our locomotor task period (i.e., $1∕150\mathrm{s}$ ). This large peak would express the changes in skin blood flow as well as brain activation depending on the task being performed. Some peaks were found at around $0.08\mathrm{Hz}$ [low frequency oscillations, (LFO)] and $0.04\mathrm{Hz}$ [very low frequency oscillation, (VLFO)]. LFO is probably related to blood pressure variations, heat-rate oscillations, and Mayers waves. (c) Power spectrum (0.2 to $2\mathrm{Hz}$ ) of Fig. 4a calculated by Fourier transform. The cardiac pulsation frequency in our locomotor task was found to be in the frequency range from 1.5 to $1.8\mathrm{Hz}$ , but the respiration frequency was not found to be in the frequency range from 0.2 to $0.5\mathrm{Hz}$ .

Fig. 5

Relevant color maps of a mixing matrix corresponding to independent components of $\mathrm{\Delta }\mathrm{oxy}$ [defined by Eq. 4] separated by MS-ICA. The component number does not have any special meaning. Each color map represents the activated spatial distribution of the corresponding component, and the distribution corresponds to the arrangement of fNIRS channels shown in Figs. 1c and 1d. The broken lines show sections of each channel. Note that these maps are spatial patterns without temporal information. The activated spatial distributions of components 1 and 3 were localized in the areas corresponding to channels 30 and 36, into which white noise was input. The activated spatial distributions of components 2, 4, 6, 7, 9, and 11 show some characteristic patterns. Interestingly, the activated spatial distribution of component 2 revealed an overall increase in the channels.

Fig. 6

Relevant independent component signal of $\mathrm{\Delta }\mathrm{oxy}$ [defined by Eq. 4] separated by MS-ICA. The component number corresponds to that of Fig. 5. Note that these graphs are temporal patterns without spatial information. The results of spectral analysis for components 1 and 3 did not give any peaks (data not shown). Their components of white noise may have been separated by ICA. The results of spectral analysis for the components 2, 4, 6, 7, 9, and 11 had one or two peaks around $0.01\mathrm{Hz}$ (data not shown). Some of these components seem to be task dependent. They may be interpreted as some independent components by the difference in brain function in the cerebral blood signal.

3.3.

New Statistical Value to Find Skin Blood Flow

Interestingly, the activated spatial distribution of component 2 in Fig. 5 revealed an overall increase in the channels. We speculate that the spatial distribution of the skin blood flow signals is not localized from the viewpoint of the zone of autonomic innervation, in contrast to brain activity signals in the motor area, because the skin blood flow is easily affected by activation of the autonomic nervous system.^{32, 33, 34} Therefore, it seems that component 2 is related to changes in the skin blood flow. To evaluate the activated spatial distributions of all components, we have defined a new statistical value, the coefficient of spatial uniformity (CSU):

Figure 7 shows the trend of the CSU values of the activated spatial distributions of all components in Fig. 5 [the column data of mixing matrix $\mathbf{A}$ was calculated using Eq. 7]. The CSU value corresponding to component 2 was a highly discriminating value (2.82) compared with the values of the other components (the second highest value: 1.41). In another subject, the highest value of the CSU was also a highly discriminating value (2.74) compared with the values of the other components (the second highest value was 1.37). These results suggest that the component related to the skin blood flow can be selected through CSU values using the column data of the mixing matrix obtained by ICA.

Fig. 7

Bar graph in the coefficient of spatial uniformity (CSU) values of the activated spatial distributions of all components in Fig. 5 [the column data of mixing matrix $\mathbf{A}$ was calculated using Eq. 7]. The CSU value corresponding to component 2 was a highly discriminating value (2.82) compared with the values of the other components (the second highest value: 1.41).

3.4.

Comparison of Skin Blood Flow Signal and Independent Component Analysis Components

The independent component signal of component 2 was compared with changes in the skin blood flow measured with a laser Doppler tissue blood flow meter (Fig. 8 ). Data have been normalized to the maximum change. The wave shape and dip timing of component 2 [Fig. 8b] were very similar to those of the skin blood flow [Fig. 8a]. This was surprising, because component 2 was derived by only ICA (statistical independence), without using any information such as task timing. The time resolution of the skin blood flow signal measured with the laser Doppler tissue blood flow meter and that of the fNIRS signals was 0.20 and $0.19\mathrm{s}$ , respectively. The coefficient of correlation was 0.724, indicating strong correlation between component 2 and changes in the skin blood flow. For another subject, the coefficient of correlation was 0.789.

Fig. 8

Comparison of skin blood flow signal and component-2 signal separated by ICA. Data have been normalized to the maximum change. (a) Changes in the skin blood flow measured by a laser Doppler tissue blood flow meter (corresponding to that in Fig. 3). (b) Independent component signal of component 2 separated by MS-ICA (corresponding to component 2 in Fig. 6).

3.5.

Removal of the Skin Blood Flow Signal from Functional Near-Infrared Spectroscopic Data

Time course changes in the fNIRS signals and the mapping image can be reconstructed with a particular component eliminated from the original fNIRS time course $\left(X\right)$ by substituting zero in the corresponding column of $\mathbf{A}$ in Eq. 5. The original map of the change in oxygenated hemoglobin and its reconstructed map at approximately $380\mathrm{s}$ are compared in Fig. 9. With component 2 eliminated, the activation of the premotor and medial somatosensory areas can be seen more clearly [Fig. 9b]. This suggests that ICA is useful for removing the skin blood flow artifact. Furthermore, eliminating components 1, 2, 3, 5, 8, 10, and 12 to 36 reduced the white noise localized in the areas corresponding to channels 30 and 36 [Fig. 9c]. Thus, ICA seems to be a powerful tool for improving the specificity of fNIRS mapping. However, further studies are needed to establish criteria for removing several signal components other than cerebral blood activity.

Fig. 9

Comparison of original and reconstructed maps of changes in $\mathrm{\Delta }\mathrm{oxy}$ [defined by Eq. 4] at approximately $380\mathrm{s}$ . (a) Original map of oxygenated hemoglobin. (b) Reconstructed map from which component 2, representing the component related to skin blood flow, was eliminated. The brain activity areas around the premotor and medial somatosensory were more localized after component 2 was removed. (c) Reconstructed map in which components 1, 2, 3, 5, 8, 10, and 12 to 36 were eliminated. The broken lines show sections of each channel.