1.

Introduction

Independent component analysis^{1, 2, 3, 4, 5} (ICA) is a useful method for analyzing functional magnetic resonance imaging (fMRI) or optical topography^{6, 7, 8} (OT) signals because they both contain components due to system noise and due to physiological changes other than neuronal activation-induced blood oxygenation changes.^{9, 10, 11, 12, 13, 14} ICA, however, still has some problems when used in practical applications. One of these problems is that the separated ICs need to be assigned by using another procedure for each case. In fMRI studies, we can extract neuronal activity-related ICs by calculating the correlation between the ICs and well-established model functions for neuronal activation signals. However, in OT studies, we have not obtained such model functions yet. We have developed an analysis procedure for extracting ICs related to neuronal activation separated from OT signals that does not need any model functions.

Functional brain imaging studies involve investigating the response of localized cortical activations to experimental tasks. Regional cerebral hemodynamic changes are thought to be related to localized cortical activations. However, these changes are sometimes affected by systemic hemodynamics. Any systemic hemodynamic changes occurring during experimental tasks interfere with the measured cerebral hemodynamic changes and obscure the results even if the actual neuronal activations related to the experimental tasks are localized in certain cortical areas.¹⁵ These systemic hemodynamic changes are considered global noise components when the measurement is intended to examine a localized brain function. When the amplitudes of these components are negligible, the interference caused by them may not be a serious problem. However, in experiments using motor tasks that affect systemic hemodynamics, the measured hemodynamic changes may contain systemic changes, such as noise that are not negligible. Nonnegligible systemic changes may also occur in experiments involving cognitive tasks that may affect the autonomic nervous system. Any analysis of these noisy OT signals must take this effect into account and treat the noise as appropriately as possible. Using a simple bandpass filter is not the best way to remove signals caused by systemic hemodynamic changes because the signals always fluctuate across a wide range of frequencies, including fluctuations due to hemoglobin changes caused by cortical activation.^{12, 13}

There are three stages in the analysis procedure we propose (Fig. 1 ). In the first stage, we use ICA to separate the signals measured by OT, which are a combination of source signals. In the second stage, ICs are selected based on mean intertrial cross-correlation (MITC), which is a measure of appropriateness for task-related ICs (TR-ICs). If a given IC is reproducible over repeated trials, it can be regarded as a TR-IC. However, both systemic changes and cortical activation changes are sometimes TR changes.^{10, 11, 13} Therefore, in the third stage, we first categorize TR-ICs into two clusters using a $k$ -means clustering method. One of the two clusters is expected to consist of TR neuronal activation-related ICs (TR-AICs) and the other of TR noise ICs (TR-NICs). This method works well when the waveforms of the two types of TR-ICs are significantly different. Then, after categorizing the TR-ICs into two clusters, we classify the clusters. We use locality of brain function to determine whether a cluster corresponds to TR-AIC or TR-NIC. If the cluster indicates an explicit locality, then it can be classified as TR-AIC. Such an approach will be also be useful for model-based analysis methods because it will provide helpful information for constructing empirical models for task-induced cerebral hemodynamic changes.

Fig. 1

Schematic diagram of data analysis. Participant data was separately analyzed in first and second stages. Then, TR-ICs, results from second stage, were collected from all participants and classified in third stage.

There have been numerous fMRI studies that used the ICA method to separate mixed activation signals (i.e., signals containing activation signals and noises). ^{16, 17, 18, 19, 20, 21, 22, 23, 24, 25} For example, Bartels and Zeki¹⁶ reported using a spatial-ICA method. To clarify whether ICs were appropriate, they checked the spatial and temporal attributes of ICs. They also determined the relationships between ICs using the IC’s temporal aspects.

In OT studies, new analysis approaches that treat noises in observed signals have been developed.^{15, 26, 27, 28, 29} Shroeter ²⁹ and Plichta ²⁸ used a model-based analysis approach based on a general linear model to extract cortical activation changes from noisy OT signals. Zhang ¹⁵ proposed using an analysis procedure without using a hemodynamic model. They assumed that the signals measured in the resting period contain only systemic or other noises and filtered out noises in the activation period by eigenvector-based filtering to reduce physiological interference in diffuse optical imaging. Morren ²⁷ demonstrated the potential of the ICA method in extracting fast neuronal signals in the motor cortex. For the data from finger tapping at $\sim 2.5\text{times}$ the heart beats, they used a second-order blind-identification algorithm, which is the same as the ICA algorithm we use. They found that the power of selected ICs at a given finger-tapping frequency was larger in a task period than in a rest period in 9 out of 14 subjects. Akgul ²⁶ also applied the ICA method (FastICA³⁰) to extract cognitive activity-related waveforms. They checked the correlation coefficient (CC) between separated ICs and a Gamma function model for brain hemodynamic response (BHR). The CC for an oddball task period was larger than the CC for a rest period. They concluded that it was plausible that waveforms extracted using the ICA method were related to BHR. In our study, to demonstrate the general applicability of the ICA method, we analyzed the data for a block-design paradigm and evaluated the results without using the hypothetical hemodynamic model function.

2.

Method

2.2.

Data Acquisition

All data were measured using an ETG-100 (Hitachi Medical Corporation, Japan), which had 24 measurement positions. The wavelengths of the irradiation light were 780 and $830\mathrm{nm}$ . The light was modulated at different frequencies. The signals for all measurement positions were sampled simultaneously at a frequency of $10\mathrm{Hz}$ . Twelve probes used to measure responses were arranged on a 3-by-3 probe holder (Fig. 2 ), with one probe holder on each side of the participant’s head. The center of the probe holder was C3 or C4 (as defined by the international 10–20 system³²), which is where the sensorimotor cortex is located.

Fig. 2

Probe configurations over left and right sensorimotor areas centered on locations C3 and C4, respectively.

2.4.

Analysis Procedure

2.4.1.

Signal preprocessing

We analyzed all the signals using a software we developed called the Platform for Optical Topography Analysis Tools, a plug-in-based analysis platform that runs on Matlab® 6.5 (Math Works, Inc.). We preprocessed all the data in the following way. First, we passed the data through a bandpass filter with a bandwidth of 0.008 to $0.5\mathrm{Hz}$ with a rectangular window. Second, we detected motion artifacts based on the criteria that the signal change between two successive samplings is $>0.3\mathrm{mM}\mathrm{mm}$ or that the support of the signal values for each measurement position is $>1.0\mathrm{mM}\mathrm{mm}$ . The measured positions of a motion artifact were excluded from the analysis even if the artifact was found only once.

2.4.2.

Signal separation

We conducted principal component analysis (PCA) on the obtained data to denoise the data set and to truncate the data dimension as the first step in signal separation. The original data were collected from 24 measurement positions that took 3000 samplings (at a sampling frequency of $10\mathrm{Hz}$ ). The PCA truncation neglects components that contribute $<1.0\%$ of the total energy. Thus, the component numbers of all the data sets were reduced from 24 to $\sim 18$ .

In the second step, we applied TDD to the PCA denoised data set separately for each participant. We set the time delay between 0 and $15\mathrm{s}$ in $1.5\text{-}\mathrm{s}$ steps. In the third step, we selected TR-ICs from all the ICs from all participants. The criterion for detecting TR-ICs was intertrial reproducibility, which was measured by using the MITC obtained from the time course data for ICs alone. This was not a model-based method, but rather a data-based one. The MITC is the average value of the interblock correlation coefficients over all the pairs of blocks. In this study, the number of pairs is 10 because the number of blocks is 5 for each session $({}_{5}C_2=10)$ . To select TR-ICs from all the ICs, in our group analysis, we set a criterion that MITC was $>0.2$ . This threshold was determined for each experimental condition based on the distribution of the MITC values. As shown in Fig. 3 , for all the experimental conditions the MITC value symmetrically distributes around zero within a range of $-0.2$ to 0.2, where the corresponding IC does not seem to be task related. We normalized all TR-ICs (mean 0, $\mathrm{s.d.}=1$ ) for further group analysis. Note that we applied smoothing to all trial data with a Gaussian kernel in order to calculate MITC (FWHM $5\mathrm{s}$ ).

Fig. 3

Histograms of MITC distributions for LT, RT, and NT. (a), (c), and (e) represent oxy-Hb, and (b), (d), and (f) represent deoxy-Hb for LT, RT, and NT, respectively.

Another criterion for finding TR-ICs is a model-based one that checks cross-correlation between the ICs and the model function. This method relies entirely on the accuracy of the model. The problem is that such accuracy has not yet been established for OT signals. However, OT signals are clear enough compared to fMRI signals for their wave-form characteristics to be determined. Therefore, it is possible to select TR-ICs based on the reproducibility of waveforms for repetitive trials. Such reproducibility is a basic concept in functional brain imaging. This criterion is more robust than the model-based one.

2.4.3.

Evaluation of task-related independent components

Assuming that TR-ICs were of two types, we classified them into TR-AIC and TR-NIC clusters. We also assumed that time courses of two types of TR-ICs differed and used a $k$ -means clustering method. The distance measure in the clustering process was based on the Pearson correlation coefficient. To classify signals using waveforms rather than amplitudes, we normalized the components and determined the sign of the IC so that the correlation between the boxcar waveform and the IC became positive. The data for the three experimental conditions [left-handed finger tapping (LT), right-handed finger tapping (RT), and no-task (NT)] were analyzed separately.

To assign classified clusters, we examined the “spatial” attribute of each IC derived from temporal ICA, which is based on only the time structure of signals and not on the spatial attribute. The activation area is thought to be localized in a certain cortical subdivision; therefore, if a TR-IC is definitely localized, then the TR-IC can be classified as a TR-AIC. In this study, however, because we focused only on the sensorimotor cortex area, the measurement area was not large enough to determine the general locality. Therefore, we adopted a laterality index as a measure of locality. If the spatial attribute of a TR-IC represents the definite contralaterality in a one-hand finger-tapping task, the TR-IC belongs to a TR-AIC. This criterion was adopted for determining which of the two clusters corresponds to the set of TR-AICs. We used the weighted laterality index (wLI)^{37, 38} to quantify the laterality of each TR-IC. The wLI indicates the degree of laterality from $-1$ (completely lateralized to the right hemisphere) to 1 (completely lateralized to the left hemisphere). To calculate the wLI, we obtained the spatial weight $W^{\prime }$ of the TR-IC under consideration for each measurement position by setting $W^{\prime }$ to zero if the spatial weight $W$ was smaller than half the value of the maximum of $W$ . Using $W^{\prime }$ , we defined wLI as

Eq. 1

$$\mathrm{wLI}=\frac{\sum W_{\mathrm{L}}^{\prime }-\sum W_{\mathrm{R}}^{\prime }}{\sum W_{\mathrm{L}}^{\prime }+\sum W_{\mathrm{R}}^{\prime }},$$

where $W_{\mathrm{L}}^{\prime }$ and $W_{\mathrm{R}}^{\prime }$ are two parts of $W^{\prime }$ , whose measurement positions are in the left and right hemispheres, respectively. The sums were taken over measurement positions in each hemisphere.

3.

Results

Figure 4 shows oxy hemoglobin (oxy-Hb) changes for LT. The time courses for clusters 1 and 2 continuous and block averaged over five trials are shown in Figs. 4a, 4c, 4b, 4d, respectively. Clear task-related changes can be seen with a few-seconds delay. The same tendency is found in the results for RT (Fig. 5 ).

Fig. 4

Time courses in clusters for LT for oxy-Hb. All ICs in each cluster are shown as gray dotted lines. Mean of ICs in cluster is shown as solid black line. In block-averaged data [(b), (c)], thin solid black lines indicate standard deviation of ICs in cluster. The task period is shown as a gray rectangle. Amplitudes of all ICs are normalized: (a) Continuous time courses in cluster 1, (b) block-averaged time courses in cluster 1, (c) continuous time courses in cluster 2, (d) block-averaged time courses in cluster 2.

Fig. 5

Time courses of clusters for RT for oxy-Hb.

In cluster 1 [Fig. 4a and 4b], there is a small overshoot of a few seconds after the onset of each task and a few seconds after it ends, which is often found in the results of motor task experiments. Steady-state periods appear in the middle of the trial but are not completely stable. The signals increase continuously and converge at a peak a few seconds after the end of the task. After the peak, the signals fall rapidly to the baseline.

Cluster 2 also shows task-related changes. However, it shows no steady-state period, and the time structure is an N-shaped pattern with its inflection points appearing a few seconds after the start and end of the task period.

The amplitude power ratios of TR-ICs are shown in Table 1 . The mean amplitude power ratios, which are averaged values of power ratios for all participants for each experimental condition, are shown in Fig. 6 . For evaluation of the power contribution ratio of selected TR-ICs to all separated ICs, we defined the amplitude power ratio with the following equation:

Fig. 6

Averaged power ratios for oxy-Hb changes. Error bar indicates standard deviation over ICs in cluster. Averaged power ratio for LT-cluster 1 is larger than that for RT-cluster 1 $(p<0.015)$ .

Table 1

Amplitude power ratio of TR-ICs for all 30 participants in percent.

Table 2

Derived values averaged for ICs in cluster.

To assign the two clusters, we checked the spatial attributes of the clusters by introducing wLI. The wLI values are $-0.47\pm 0.64$ and $-0.11\pm 0.74$ for LT-clusters 1 and 2, respectively (Table 2). The wLI of LT-cluster 1 shows clear contralaterality ( $p<0.01$ ; hypothesis is $\mathrm{wLI}=0$ ), but cluster 2 does not $(p=0.55)$ . The same results are also found in the wLI values for RT clusters.

Results for deoxy hemoglobin (deoxy-Hb) are shown in Figs. 7 and 8 and Tables 1, 2. The extracted time courses are almost the same as those for oxy-Hb (Figs. 7 and 8). The amplitude power ratios of deoxy-Hb differ significantly between clusters 1 and 2 for deoxy-Hb but not for oxy-Hb (Table 1). The averaged MITCs for deoxy-Hb are similar to those for oxy-Hb (Table 2). The values of wLI show contralaterality for both LT and RT that is statistically significant in LT $(p=0.01)$ but not in RT $(p=0.20)$ .

Fig. 7

Time courses of clusters for RT for deoxy-Hb.

Fig. 8

Time courses of clusters for LT for deoxy-Hb.

The mean-amplitude-power ratios of TR-ICs are $3.5\pm 7.2\%$ and $1.3\pm 2.6\%$ for oxy- and deoxy-Hb concentration changes, respectively, for NT (Table 1). These ratios are smaller than those for LT and RT. The numbers of TR-ICs found in all extracted ICs from all participants are also smaller than those in LT and RT. For oxy-Hb concentration changes, the numbers of the TR-ICs found are 55, 51, and 14, while the total number of ICs are 522, 503, and 476 for LT, RT, and NT, respectively. For deoxy-Hb concentration changes, the numbers of the TR-ICs found are 44, 42, and 8, while the total numbers of ICs are 540, 512, and 512 for LT, RT, and NT, respectively. Figure 9 shows the time courses of TR-ICs for oxy-Hb concentration changes for NT.

Fig. 9

Time courses of clusters for NT for oxy-Hb.

4.

Discussion

5.

Conclusion

We proposed an analysis procedure for extracting ICs related to neuronal activation separated from OT signals without using any model functions. We used a TDD for signal separation, a MITC as a criterion for selecting TR-ICs, and a clustering method for classifying TR-ICs into TR-AICs and TR-NICs, such as physiological changes. The data set we used was obtained from one-handed finger-tapping experiments. We found that a cluster of TR-AICs showed distinct contralaterality, whereas another cluster, thought to be of TR-NICs, did not. We obtained TR-AICs without using model-based analysis. We conclude that the analysis procedure we proposed in this study is useful for analyzing OT data.