Detection of fast neuronal signals in the motor cortex from functional near infrared spectroscopy measurements using independent component analysis

Fast changes, in the range of milliseconds, in the optical properties of cerebral tissue are associated with brain activity and can be detected using noninvasive near-infrared spectroscopy (NIRS). These changes are assumed to be caused by changes in the light scattering properties of the neuronal tissue. The aim of this study was to develop highly sensitive data analysi algorithms to detect this fast signal, which is small compared with other physiological signals. A frequency-domain tissue oximeter, whose laser diodes were intensity modulated at 110 MHz, was used. The amplitude, mean intensity and phase of the modulated optical signal were measured at a sample rate of 96 Hz. The probe, consisting of four crossed source detector pairs was placed above the motor cortex, contralateral to the hand performing a tapping exercise consisting of alternating rest and tapping periods of 20 s each. An adaptive filter was used to remove the arterial pulsatility from the optical signals. Independent component analysis allowed further separation of a signal component containing the fast signal. In nine out of 14 subjects, a significant fast neuronal signal related to the finger tapping was found in the intensity signals. In the phase signals, indications of the fast signal were found in only two subjects.

approach to the study of neurovascular coupling (VILLRINGER and CHANCE, 1997), combining sub-centimetre spatial localisation with millisecond temporal resolution ). An important advantage of NIRS is that it is non-invasive and can be used at the bedside. Owing to the limited penetration of light into living tissue, however, only superficial layers of the brain can be reached with NIRS .
The present paper deals with the detection of the fast neuronal signal in the motor cortex. As a finger-tapping exercise was used as stimulation, the fast signal is correlated with the tapping frequency, it is, however, very small compared with other physiological signals, such as arterial pulsatility (due to systole and diastole). Therefore highly sensitive data analysis algorithms have been developed, in a previous study, it has already been shown that the fast signal can be detected in the motor cortex using frequency-domain NIRS . The detection method was based on the cross-correlation function of the optical data with the tapping signal (finger movement) over a measurement period of 5 min.
The main contribution of the present paper is the use of more advanced signal processing techniques allowing the separation of a signal component containing the fast signal. Using adaptive filtering and independent component analysis (ICA), signal components are extracted in which the individual stimulation periods of 20 s can be identified.

Subjects
Three female and 11 male healthy adult volunteers (between 19 and 55 years) were included in this study. All except one were right handed. Written informed consent was obtained from all subjects prior to the measurements.

Measurements
An instrumental low-noise system was developed for the NIRS measurements MORNER et al., 2002). A two-wavelength frequency-domain tissue oximeter*, whose laser diodes were intensity modulated at 110 MHz, was used. The light of four laser diodes, at either 758 nm (sources 1 and 3) or 830 nm (sources 2 and 4), was combined to provide a higher light intensity at each source location. The power output of each laser diode was 5 mW (laser class 3a). All sources and detectors were arranged on a circle with a diameter of 3.2 cm, as shown in Fig. 1. So that the light of the four source locations could be distinguished, the laser diodes were multiplexed. At each of the four detectors (A,B,C,D), the amplitude (AC), mean intensity (DC) and phase (PH) of the modulated optical signal were measured. The sample rate for a complete cycle, which measures all 48 parameters (AC, DC and PH for 16 sourcedetector pairs), was 96.1539Hz. Representative fragments of DC and PH signals are shown in Figs 2a and b. it can already be seen that the PH signals are much more noisy than the DC signals.
in addition to these optical signals, the arterial pulsations, the respiration and a signal synchronous with the finger tapping were measured. The arterial pulsations were acquired by means of a pulse oximeter t attached to the left hand, whilst the tapping was performed with the right hand and vice versa. The respiratory signal was obtained from a breathing belt fixed around the chest. Using a strain gauge ~, a signal proportional to the tension in the belt was derived. An electronic device

Protocol
The sensor, consisting of the four crossed source detector pairs, was placed above the motor cortex (C3 position according to the international 10/20 system (JASPER, 1957)), contralateral to the hand performing the finger-tapping exercise. The tapping frequency was set, with a metronome, at approximately 2.5 times the heart rate of the subject to avoid the influence of harmonics of the arterial pulsations on the fast neuronal signal. The measurement protocol consisted of a 1 min pre-exercise period (baseline), a 5 min stimulation run with alternating rest and tapping (with the whole hand) periods of 20 s each (corresponding to a frequency of 0.025 Hz), and a 1 min after-exercise period (baseline). For two subjects, the stimulation runs had different durations. Overall, there were 14 subjects: 12 with stimulation runs of eight tapping periods, one with 16 tapping periods and one with six tapping periods. The subjects were supine in a quiet, dark room and wore safety goggles during the whole measurement. The protocol was approved by the institutional Review Board of the University of illinois at Urbana-Champaign ORB 94125).

Signal processing
The fast neuronal signal is so small compared with other physiological signals that it cannot be directly observed in the NIRS signals, in fact, the largest alterations in the optical signals are caused by the arterial pulsatility due to systole and diastole. This pulsatility was first reduced in each NIRS signal using an adaptive filter. However, even after filtering, the fast signal can hardly be detected in individual NIRS signals. Apparently, a multichannel technique is needed to exploit the information in different NIRS signals (light bundles) simultaneously. Therefore independent component analysis (ICA) was applied to the set of all filtered NIRS signals (per subject) to extract a signal component where the fast signal could be detected more clearly.

Adaptive filtering
The physiological noise due to the arterial pulsations was removed from the NIRS signals using an adaptive filter. Two different filters were compared. First, a standard interference canceller was used, as shown in Fig. 3, where the adaptation is based on the normalised leastmean-square (NLMS) algorithm (HAYKIN, 1996). The primary signal, which serves as the desired response for the adaptive filter, is, in our case, the NIRS signal containing the signal of interest (tapping) and noise (arterial pulsations). The pulse oximeter signal, which contains the arterial pulsations but not the tapping signal, is used as input to the adaptive filter (noise reference) to reduce the noise in the NIRS signal. An adaptive filter with 300 taps was used, and the adaptation constant of the NLMS algorithm was set to 0.1.
in some measurements, however, the pulse oximeter signal was saturated, resulting in a distorted waveform, resembling a square wave, of the arterial pulsations, in those cases, a filter was used that only uses the NIRS signals and not the pulse oximeter signal. This filter is described in detail elsewhere (GRATTON and CORBALLIS, 1995). in short, it extracts a mean shape of the pulse by screening each trace separately for pulses, whose period is adjusted before the averaging. This mean shape corresponds to the best estimate of each pulse and is used to remove each pulse from the data. This again requires adjusting the period of the mean shape to that of each pulse and scaling the shape by linear regression, in the remainder of this paper, this filter is abbreviated to HF.

Independent component analysis (ICA)
ICA is a technique that allows blind recovery (ICA is also called 'blind source separation') of statistically independent primary. signal The following linear statistical model is assumed (HYV,~RINEN et al., 2001;COMON, 1994): where x and n are vectors of size m (the number of light bundles) representing the observations and the disturbances (in the form of additive noise), and s is a vector of size n representing the n sources. The goal is to identify the m x n mixing matrix A and/or the source vector s, only from the corresponding observations x. The basic assumption used to solve this problem is that the components of s are mutually statistically independent, as well as independent of the noise components. Assuming that the tapping is statistically independent of the other physiological sources (heart rate, respiration), ICA should be able to separate them. The ICA problem is usually solved by a two-stage algorithm. in the first step, an n x m matrix W is calculated such that the transformed observations z = Wx are 'white', i.e. its components are uncorrelate,d, and their variances equal unity (hence, E[zz*] = l, where z stands for the hermitian transpose of z). The whitening matrix W can be determined using only second-..... order statistics (i.e. from the covarlance matrix R x = E[xx ] of the observations). This prewhitening step, which amounts to a principal component analysis (PCA), can also be used to reduce the dimensions of the data.
in a second step, a unitary n x n matrix U is calculated such that A = pinv(W)U, where piny(W) denotes the Moore-Penrose pseudo-inverse, in 'general' ICA algorithms, the matrix U is determined from higher-order statistics, (COMON, 1994;CARDOSO, 1999;DE LATHAUWER, 1997). For temporally correlated signals, however, the time dependence structure can be exploited, and U can also be retrieved from covariance matrices at non-zero lags (hence, using only second-order statistics) (BELOUCHRANI et al., 1997;ZlEHE and MOLLER, 1998;MOLGEDEY and SCHUSTER, 1994).
As the optical signals are clearly temporally correlated, we used an algorithm that exploits this structure, namely SOBI (BELOUCHRANI et al., 1997). The SOBI algorithm is based on the simultaneous diagonalisation of a set of covariance matrices at different lags. For this study, we used ten covariance matrices at lags of 10, 20,..., 100 samples. Experiments with various other lags showed that these values could yield the best results for our data. The performance of SOBI was compared with that of several other algorithms (JADE (CARDOSO, 1999); fastlCA (HYV,~RINEN and OJA, 2000); infomax (BELL and SEJNOWSKI, 1995); TDSEP (ZlEHE and MtiLLge, 1998)), including some 'general' ICA algorithms, for this application, but these did not give any better results.

Results
From the theory of photon migration, it is known that DC and AC signals contain the same information. Therefore, as the signal-to-noise ratio (SNR) of the AC signals is smaller than that of the DC signals , only the DC and PH signals were analysed in this study. Fig. 2a shows the CDC3 signal (DC, detector C, source 3) in subject 3 (according to Table 1). This signal is representative of the other DC signals, although the signal strength depends on the distance between source and detector and differs between subjects. The dominant frequency is due to the arterial pulsations. As shown in Fig. 2b, the PH signals are much more noisy.
The power spectra of representative DC and PH signals (ADC4 and APH4) were estimated using the Welch method (KAY, 1988), with non-overlapping windows of 15 000 samples (156s), and are shown in Fig. 4. A signal of the longer measurement of subject 9 was used for Fig. 4 as this allowed us to average several periodograms of such long windows.
in the power spectrum of ADC4, the dominant frequencies between 1 Hz and 5 Hz are the heart rate (arterial pulsations), around 1.4 Hz, and its harmonics. As the heart rate is not perfectly constant during the whole measurement, the peaks are smeared out. Atthe tapping frequency (3.55 Hz), no peak can be distinguished without prior knowledge about the tapping frequency. There are also peaks in the lower frequency range (< 1 Hz). The peak around 0.3 Hz is due to the respiration. This peak is, however, not always present in the power spectra of the optical signals. The peaks around 0.12 Hz are related to Mayer waves, which have been observed by many researchers, although their origin has not yet been exactly determined (HUDETZ et al., 1998). The peak at 0.025 Hz, corresponding to a period of 40 s, is due to the slow haemodynamic signal. This is a first indication that the measurements were taken over the activated region of the brain. In the power spectrum of APH4, only a peak at the heart rate can be observed, which is usually smaller than the corresponding peak in the DC signals. In many cases, the noise in the phase signal is so large that even this peak cannot be detected.

Reduction in arterial pulsatility using adaptive filtering
The optical signals were first filtered with an eighth-order bandpass Butterworth filter, with a passband between 1 and 10 Hz. in this way, the slow fluctuations due to respiration and Mayer waves were removed from the optical signals. Consequently, the arterial pulsations were removed using an adaptive filter, as described in Section 2.4.1. Comparing both adaptive filters (NLMS and HF), it was noted that the NLMS filter yielded better results (in terms of reduction of the peaks at the heart rate, as observed in the power spectrum), except in some cases where the pulse oximeter signal was heavily saturated. For the latter signals, the HF filter from GI{ATTON and CORBALLIS (1995) was used. The CDC3 signal of subject 3 in Fig. 5a after filtering (bandpass and adaptive NLMS) is shown in Fig. 5b. Fig. 5a also shows the CDC3 signal after adaptive filtering only (no bandpass filtering) to illustrate clearly the reduction in arterial pulsations.
The effect of the adaptive filtering is better visible when the power spectra of both signals are compared, as shown in Fig. 6. Some residual of the arterial pulsations can be observed in the power spectrum of the filtered CDC3 signal as peaks at the second and third harmonics of the heart rate (around 2.6 Hz and 3.9 Hz respectively). These power spectra were estimated using the Welch method with non-overlapping windows of 1920 samples (20 s), corresponding to the tapping and rest periods. The peaks at the heart rate and its harmonics are almost completely removed in the filtered signal. The peak at the tapping frequency, however, is still so small that it would be difficult to distinguish it without prior knowledge about the tapping frequency.  To determine whether the observed peaks in the power spectra were actually caused by the stimulation, the power at the tapping frequency during consecutive tapping and rest periods was compared. The power spectra were calculated in this case for each tapping and rest period, with three overlapping windows, whose length was half the considered period. For each stimulation period, the tapping period was determined from the spectrum of the tapping (movement) signal; the same frequency was used to calculate the power during the following rest period. The Wilcoxon rank-sum test was used to determine whether the power at the tapping frequency during the tapping periods was significantly different from that during rest periods (LEHMANN, 1975). This test, which calculates the probability that the given data (or more extreme) would be observed if both populations had the same mean, allows us to compare and summarise the results of the different measurements in a simple way.
The results are summarised in Table 1, which gives the lowest p-value among all channels for each measurement, together with the number of channels with ap-value less than 0.05. in some measurements, one or two optical signals were saturated and therefore excluded from the analysis. For each significant signal, it was verified that the difference was in the right direction, i.e. that the power during the tapping periods was larger than that during the rest periods, in some measurements, where the peaks at the 2nd or 3rd harmonic were more clear in the power spectrum than the peak at the tapping frequency itself, the respective harmonic was used for the detection, as indicated in Table 1.
In 49 DC signals out of 206 (24%), the power at the tapping frequency during the tapping periods was significantly different from that during the rest periods, at a level ofp < 0.05. In eight subjects (out of 14), more than one channel withp<0.05 was found. For the smallest p-value among all channels for each measurement, the significance level was lowered to p <0.003, according to multiple testing corrections. Note that this is a conservative approach based on Bonferroni inequalities, although the channels are not independent, as they interrogate the same tissue, and that it is actually not expected that a fast signal would be seen in all channels. In nine (out of 14) subjects, a significant DC channel was found at thep < 0.003 level; seven of these correspond to subjects where more than one channel withp < 0.05 was found.  Fig. 6 Power spectra of signals" in Fig. 5. (-.-.) Spectrum of unfiltered CDC3 signal is" dominated by peaks" at heart rate and its" harmonics, owing to arterial pulsations. These peaks" are almost completely removed in (---) power spectrum of filtered (bandpass and adaptive NLMS) CDC3 signal, but peak at tapping frequency is" still very small. Dominant peaks" in power spectrum of second independent component IC2 (calculated from set of all filtered DC signals" of subject 3) are at tapping frequency (3.35Hz) and its" second harmonic (6.7Hz). Note: spectra were shifted vertically to produce clear figure In ten PH signals out of 206 (5%), ap-value below 0.05 was observed. One channel (subject 2) was significant at the level of p < 0.003. The two channels in this subject that were significant at thep < 0.05 level were also significant at thep < 0.05 level in the corresponding DC channels (same source-detector pair), it was also noted that, in subject 9, when the 3rd harmonic was analysed (not shown in Table 1), six DC channels and eight PH channels were significant at the p < 0.05 level; five of those channels were significant in both DC and PH.

Extraction of fast neuronal signal using ICA
To separate the fast signal further, ICA was applied to the optical signals over the whole exercise epoch. The optical signals of each measurement were filtered, with either NLMS or HF, and subsequently SOBI was applied to the set of all nonsaturated optical signals of that subject, in some measurements, one or two optical signals were saturated and therefore excluded from the analysis. The number of independent components (iCs) calculated was always taken to be equal to the number of sources. Fig. 5c shows one of the iCs, calculated from the DC channels of subject 3, in which the tapping frequency is present. As shown in Fig. 6, the dominant peaks in the power spectrum of this IC are at the tapping frequency (3.35 Hz) and its second harmonic (6.7 Hz). The other iCs contain the (residual) arterial pulsations and/or noise. The fact that these peaks are caused by the finger tapping is clear from Fig. 7, which shows the power at the tapping frequency during consecutive tapping (and rest) periods.
The Wilcoxon rank-sum test was used to determine whether the power at the tapping frequency during the tapping periods was significantly different from that during the rest periods. For the IC shown in Fig. 5c, the p-value is 0.0002. In Table 2, the lowest p-value among the iCs is given for each measurement. Considering the DC signals, in nine subjects a significant IC was found at thep < 0.003 level, it is noted that, for all these subjects, the significance of the iCs is higher than that of the DC channels (lower p-value), reflecting the improved separation of the tapping signal due to ICA. Concerning the PH signals, in four subjects the lowestp-value among all iCs was lower than 0.05, including the subject with a significant PH channel at the p<0.003 level. Considering the 3rd harmonic of subject 9, which has eight PH channels with a p-value below 0.05, the lowest p-value observed among the iCs was p = 0.0001.
The independent component containing the fast signal also allows the tapping periods to be identified individually. Fig. 8a shows the spectrogram of the tapping signal (finger movement), and Fig. 8b shows the spectrogram of the IC containing the fast signal. The spectrograms were calculated as the short-time Fourier transforms (STFTs) over non-overlapping windows of 4 s. All tapping periods, except the first, can clearly be distinguished in the spectrogram of the IC. Even the shift of the tapping frequency in each tapping period and between tapping periods can be observed. Fig. 9 shows the power at the tapping frequency, as in Fig. 7, but for 'arbitrary' (not perfectly coinciding with tapping and rest periods) time windows of 4 s. The tapping frequency was determined from the spectrogram of the IC (optical signal) itself, by the selection, for each instant in time, of the maximum of the STFT over the frequency range of interest. To determine, in a more objective way, which tapping periods can be detected, the following simple procedure was followed. The mean plus twice the standard deviation of the power at the tapping frequency during the baseline period (the first minute of the measurement, without stimulation) was used as the threshold level. A tapping period (of 20 s) was considered as detected when the threshold was exceeded three or more times (out of five, for non-overlapping windows of 4 s) during that tapping period. As shown in Fig. 9, all except the first tapping period in the spectrogram of Fig. 8 are detected in this way. Also, for the other measurements, the tapping periods detected in this way agree very well with those that can be distinguished by eye in the corresponding spectrograms. As this procedure does not use the tapping signal in any way, it allows us to quantify which tapping periods can be detected using only the optical signals, as shown in Table 2. Considering the nine subjects with significant iCs at thep < 0.003 level (calculated from the DC channels), 67 out of 78 tapping periods (85%) have been detected individually.

Discussion
Considering the DC signals, in nine subjects an IC was found in which the power during tapping was significantly larger than during rest periods (subjects 1-5, 7, 9, 10, 12). The significance level ofp < 0.003 that we used is a conservative value, taking into account multiple testing corrections. Furthermore, the nine subjects include all eight subjects where several significant DC channels at thep < 0.05 level were found; the remaining subject with an IC with ap-value below 0.003 had one DC channel that was also significant at the p < 0.003 level, in six of these nine subjects, the p-value was actually the lowest possible value that can be obtained with the rank-sum test for the given sample size (eight or six data points), meaning that the power at the tapping frequency during any tapping period is larger than the power during all rest periods. Considered together, these findings show that the iCs considered contain a fast signal relative to the finger tapping, in these nine subjects, where the fast signal was detected in the DC signals, 67 out of 78 tapping periods could be identified individually in the spectrogram of the IC containing the fast signal. As the fast signal is expected to be very localised, the most likely reason for not finding the fast signal in all subjects is that the sensor, which is relatively small, was not precisely in the right position. There was one subject with a significant PH channel at the p < 0.003 level. The fact that the two channels in this subject that were significant at the p < 0.05 level were also significant at the p < 0.05 level in the corresponding DC channels (same sourcewere detector pair) suggests that the observed low p-values are not due to chance. The same could be assumed for the 3rd harmonic in subject 9, where five (out of 16) channels were significant in both DC and PH at thep < 0.05 level, and the lowest p-value among the iCs was p = 0.0001. in general, however, the fast signal could not be reliably detected in the PH signals, probably because there is too much noise in the PH signals. Our statistical analysis shows that the described fast signal in the DC signals is not due to random noise. Furthermore, it cannot be a movement artifact caused by direct mechanical coupling between the tapping and the optical sensors because it is localised. A previous study on a subset of the measurements analysed in the present paper showed that the localisation of the fast signal was in reasonable agreement with the slow signal . Therefore we conclude that the observed fast signal originates from neuronal activity.
Relatively few other studies describing a fast neuronal signal detected non-invasively in humans using NIRS have been published (GRATTON et al., 1995a;R[NNE et al., 1999;STEINBRrNK et al., 2000;DESOTO et al., 2001;GRATTON and FABIANI, 2003;WOLf et al., 2003). Gratton et aL report increases in the phase of a frequency-domain NIRS instrument during various types of stimulation (GRATTON et al., 1995a;RINNE et al., 1999;DESOTO et al., 2001). Another group, from Berlin, failed to detect a change in the phase in the visual cortex using a similar experimental set-up (SYR~ et al., 2000). According to Gratton, this can be explained by the fact that the study did not include any systematic mapping of occipital areas, nor contralateral control of stimulated areas, nor averaging across subjects (GRATTON and FABIANI, 2003b).
The Berlin group also conducted a study using continuouswave NIRS instrumentation, where a decrease in the intensity of the reflected light was observed following somatosensory stimulation . in the data presented here (from the motor cortex), as well as in other data (from the visual cortex) recorded in the same laboratory with frequencydomain NIRS, a significant fast signal was observed in the intensity but not in the phase . We only found indications of the fast signal in the phase in two out of 14 subjects, which is probably owing to the high noise in the phase signals. Based on studies conducted in single neurons (STEPNOSKI et al., 1991), in populations of neurons at the macroscopic level (RECTOR et al., 1997) and (invasively) in animals (MALONEK and GRINVALS, 1997), all authors of the cited literature assume that the fast signal originates from changes in light scattering. However, it cannot be excluded that the reported fast signals in human subjects are related to changes in light absorption. As we could not reliably detect the fast signal in the phase, we are not able to clarify this aspect.
An important methodological difference between the present work and the cited studies on the fast neuronal signal concerns the signal processing techniques. Most other studies use a folding average approach, which is not effective in our case because the tapping frequency is not perfectly constant, contrary to the stimulation frequency in the case of visual or somatosensory stimulation. Using adaptive filtering and ICA, however, we were able to separate a component from the optical signals in which individual stimulation periods can be identified. The advantage of the presented approach, which can also be applied in other cases such as visual stimulation, is that less measurement time is required per subject or more types of test can be carried out within the same subject. Furthermore, new effects, such as habituation, can be studied. Localisation based on ICA is another potential application.

Conclusions
Low-noise frequency-domain NIRS instrumentation, using four crossed source-detector pairs, was used to study the fast neuronal signal related to a finger tapping exercise. The arterial pulsatility, which was the main source of physiological noise, was first reduced using an adaptive filter. On the filtered signals, independent component analysis was applied to separate a signal component due to the neuronal activity. In nine out of 14 subjects, an independent component could be extracted from the intensity signals where the power at the tapping frequency during stimulation was significantly larger than during rest. In those nine subjects, 67 out of 78 tapping periods could be identified individually in a time-frequency plot of the independent component. In the phase signals, we only found indications of the fast signal in two out of 14 subjects. the Flemish government: FWO projects G0200.00, G0078.01, G0407.02, G0269;02, G0270.02 and G0240.99, by the Belgian federal government: DWTC (IUAP V-22) and by the EU: NICONET, INTERPRET, PDT-COIL, MRS/MRI signal processing (TMR). This research was supported by NIH grant CA57032.