BS3A.21.pdf
Biomedical Optics 2014 © OSA 2014
Characterization of the NIRS Hemodynamic Response Function with Independent Component Analysis Rodrigo M. Forti1, Andréa Alessio2, 1, Rickson C. Mesquita1 1
Institute of Physics “Gleb Wataghin”, 2 Faculty of Medical Sciences, University of Campinas (UNICAMP), Campinas, SP (Brazil) Rua Sérgio Buarque de Holanda, 777 - Cidade Universitária Zeferino Vaz - Barão Geraldo – Campinas - SP - Brazil Author e-mail address:
[email protected]
Abstract: We propose an algorithm with objective criteria to identify the hemodynamic response function by using independent component analysis during cerebral activation with NIRS. The algorithm was successfully validated in both real data and simulations. OCIS codes: (100.1455) Blind deconvolution; (200.4560) Optical data processing.
Introduction: Independent Component Analysis (ICA) is a blind-source separation technique that uses statistical independence to separate source components (ICs) from a linear mixed signal [1]. ICA has been previously used in functional NIRS data mostly to remove signal artifacts [2, 3] and to identify spatial maps related to brain connectivity [4]. However, the degree of freedom in choosing the number of ICs and their estimated amplitudes makes the technique difficult to fully characterize the hemodynamic response function (HRF) with ICA in functional NIRS (fNIRS) studies. In this work we aimed to investigate the feasibility of establishing objective criteria to identify the HRF in fNIRS activation experiments using ICA. We assumed physiological and physical constrains to find the IC related to the HRF, and investigated how many ICs must be kept after the whitening process. We tested our criteria on simulated and real datasets, and found that we can always find at least one IC that is strongly correlated to the HRF estimated by traditional methods (i.e., block averaging or deconvolution process). Materials and Methods: In order to test the feasibility of ICA to estimate the HRF, we employed our novel algorithm on two different datasets: a simulated fNIRS dataset and data from a real experiment involving a cognitive task. First, we simulated an fNIRS experiment with 28 source-detector separations (i.e., channels) at 3 cm each. For each channel we simulated a time-course, and we hypothesized that the temporal series were composed of up to three sources: 1) local brain activation, that was simulated as a Gaussian signal with an initial deep, and it was modulated by a block-design stimulation paradigm; 2) global physiology signal, that was simulated with a sinusoidal function with peak frequencies varying from 0.9 to 1.3 Hz, and; 3) extra-cortical signal, that was simulated with a Rayleigh distribution and it was also modulated by the stimulus [5]. The global physiology signal was added to all channels, while the local brain activation and the extra-cortical signals were added to a fraction of the channels, covering from 10% to 90% of all channels. Each signal was randomly added to a channel with normalized amplitudes. We also added Gaussian (white) noise to each channel, with amplitudes 0.25 and 0.5, respectively, in order to simulate instrumental and other potential sources of noise. All simulations were performed at least 7 times for analysis of statistics and systematic errors. Next, we analyzed the performance of our ICA algorithm on real data. We recruited 7 healthy subjects, all male, right-handed, with ages varying from 19 to 26 years old. Data were acquired with a commercial continuous-wave NIRS system (NIRScout, NIRx Medical Technologies, Glen Head, NY, EUA) with temporal resolution of 0.1 s. The geometry had channels with source-detector distances varying from less than 2 cm (short channels) up to 3 cm (long channels), in a total of 50 channels. We expect that the short channels exhibit extra-cortical contributions only, while the long channels would contain both extra-cortical and cortical contributions to the NIRS signal. For functional activation, we performed a simple verbal fluency task. The task consisted of one run with 5 stimulation blocks of 30 seconds each. In each block, the subject was required to verbalize words starting with a given letter (F, A, C, B, M). All procedures were approved by the Institutional Review Board of the University of Campinas, where the experiments were carried out. After acquisition, data were filtered between 0.01 Hz and 3 Hz, in order to reduce instrumental noises. The modified Beer-Lambert law [6] was used to convert the measured intensity into oxy- and deoxy-hemoglobin concentration changes (ΔHbO2 and ΔHb, respectively) at each channel. For the functional activation analysis and estimation of the HRF, we first analyzed the data with traditional block-averaging and deconvolution methods. The results were later compared to our proposed ICA approach.
BS3A.21.pdf
Biomedical Optics 2014 © OSA 2014
For ICA analysis, we performed ICA decomposition on the HbO2 and Hb time-courses with the FastICA algorithm [1]. In order to identify the IC related to the HRF, we established two physiological criteria: 1) the ratio 𝑅𝑓𝑟𝑎𝑐 of the peaks between the physiological noise region (0.8 – 1.7 Hz) and the brain activation region (0.01 – 0.25 Hz) in the frequency spectrum must be equal to or less than 0.15; and 2) the mean inter-trial cross-correlation (MITC) of the computed IC must be higher than 0.45. The MITC measures the reproducibility of the HRF across different trials/stimulus [7]. Our algorithm initially computed 2 independent components, and iteratively increased the number of ICs until the above-mentioned criteria were met. All scripts were run in Python. Results and Discussion: With the simulated dataset, the algorithm with ICA was able to find the simulated local activation (i.e., the HRF) at any given fraction of source signals in the time-series. For the case where we had no global signal in the data, the algorithm always found an IC related to the local signal (i.e., an IC that obeyed our established criteria) with only 2 independent components computed; the Pearson’s correlation coefficient between the estimated HRF and the simulated local signal was ρ = 0.99±0.2. When the local signal was mixed with a global signal, the algorithm found the HRF computing 3 ± 1 ICs. In addition, the ICA procedure was able to separate the global signal from the local signal in the estimated HRF without any prior information (Figure 1). The correlation between the estimated IC and the simulated local signal yielded ρ = 0.98±0.01. The results achieved with ICA are significantly better than the results obtained with the traditional method. The global and local signals could not be separated using traditional methods, yielding a correlation of ρ = 0.76 ± 0.08 between the HRF estimated by block-average and the simulated local signal. Our results were independent of the signal-to-noise ratio (SNR) of the data, with the ICA algorithm performing better than the traditional approach in all cases tested. Also, the algorithm correctly identified the HRF when the local activation was present in 10% or 25% of the channels, suggesting that the proposed approach is very sensitive to local perturbations, such as short finger-tapping or visual tasks. For the real experiment, we found consistent activation patterns in the temporal region of the brain across all subjects. Using the traditional approach for estimating the HRF, we found a significant increase in the amplitude of both oxy-hemoglobin and deoxy-hemoglobin concentrations around Broca’s area, in agreement with other studies of verbal fluency [8, 9]. Across all subjects, Δ𝐻𝑏𝑂2 = 0.6 ± 0.5 𝜇𝑚𝑜𝑙 and Δ𝐻𝑏 = 0.2 ± 0.2 𝜇𝑚𝑜𝑙. ICA was performed in all subjects for each concentration data independently. An example of the HRFs found for a representative subject is shown in Figure 2. For oxy-hemoglobin, our objective criteria found at least one IC that could represent the HRF. In all cases, the algorithm stopped running with Rfrac ≤ 0.14 and MITC ≥ 0.55. The IC was highly correlated to the HRF estimated with the traditional methods when only the long channels where considered (ρ = 0.8 ± 0.2), but presented a low correlation when the short channels were considered (ρ = 0.3 ± 0.7). In 2 subjects, the algorithm found two independent components that met the criteria established. After careful analysis, we noticed that one IC was highly correlated to the HRF in the long channels while the other one was similar to the HRF considering the short-channels only. That result suggests that our ICA algorithm has potential to separate cortical from extra-cortical contributions even in real human data, similarly to what we found in our simulation studies. For the deoxy-hemoglobin time-series, the ICA approach yielded lower correlation coefficients. The ICs found by the algorithm presented an average correlation of ρ = 0.6 ± 0.5 and ρ = 0.4 ± 0.6 when compared to the HRF estimated by block-average considering the long and short channels, respectively. On 3 subjects, our criteria found two independent components. In these cases, one component was highly correlated to the HRF in the short channels, while the other component did not show any significant correlations with the traditional HRF. In all cases, the algorithm stopped with Rfrac ≤ 0.08 and MITC ≥ 0.45.
Figure 1 - Example of the HRFs simulated and obtained with our criteria on the simulated data. The gray area represents the simulated stimulation period.
BS3A.21.pdf
Biomedical Optics 2014 © OSA 2014
Figure 2 - Example of the HRFs estimated with the traditional (blue and red) and our ICA (black) approach. For the traditional method, each curve represents an average over the activated region (Broca’s area) during the verbal fluency task. The gray area is the stimulation period.
Conclusions: In this work, we propose the use of ICA to find the HRF in an objective approach. By establishing two simple and reasonable assumptions, our results show that it is possible to identify at least one IC that represents the HRF, and that this IC is strongly correlated with the HRF estimated by traditional fNIRS analysis. Such correlation is more evident in the HbO2 time-series as compared to Hb curves. This result might be explained because, in our experiment, the estimated HRF did not change whether we considered extra-cortical contributions or not. In agreement with previous findings [7], our simulations showed that our algorithm could separate global systemic physiology from local brain HRF even if they were correlated with the stimulus. One limitation of the algorithm, however, is the lack of information in the amplitude of the estimated HRFs, since this parameter is arbitrary in the ICA transformation. Although this limitation makes it harder to compare different subjects, differences in the HRF dynamics are important feature of cerebral activation, and very useful for studying parametric variations, or when comparing different populations, such as disease vs. control groups, during tasks. References: [1] Aapo Hyvärinen, “Fast and Robust Fixed-Point Algorithms for Independent Component Analysis”, IEEE Trans. on Neural Networks, 10(3):626-634 (1999). [2] S. Kohno, I. Miyai, A. Seiyama, I. Oda, A. Ishikawa, S. Tsuneishi, T. Amita, K. Shimizu, “Removal of the skin blood flow artifact in functional near-infrared spectroscopic imaging data through independent component analysis”, Journal of Biomedical Optics 12(6), 062111 (2007). [3] J. Virtanen, T. Nononen, P. Meriläinen, “Comparison of principal and independent component analysis in removing extracerebral interference from near-infrared spectroscopy signals”, Journal of Biomedical Optics 14(5), 054032 (2009) [4] J. Markham, B. R. White, B. W, Zeff, J. P. Culver, “Blind Identification of Evoked Human Brain Activity With Independent Component Analysis of Optical Data”, Hum. Brain Mapp. (2009) [5] T. Takahashi, Y. Takikawa, R. Kawagoe, S. Shibuya, T. Iwano, S. Kitazawa, “Influence of skin blood flow on near-infrared spectroscopy signals measured on the forehead during a verbal fluency task”, NeuroImage 57 991-1002 (2011). [6] L Kocsis, P Herman, A Eke, “The modified Beer–Lambert law revisited”, Phys. Med. Biol. 51 N91–N98 (2006) [7] T. K. H. Sato, Y. Fuchino, T. Yoshida, H. A. M. K. A. Maki, M. Abe, N. Tanaka, “Extracting task-related activation components from optical topography measurement using independent components analysis”, Journal of Biomedical Optics 13(5), 054008 (2008). [8] T. Sutoa, M. Fukudaa, M. Itoa, T. Ueharaa, M. Mikunia “Multichannel near-infrared spectroscopy in depression and schizophrenia: cognitive brain activation study”, BIOL PSYCHIATRY 55:501–511 (2004) [9] M.Kameyama, M. Fukuda, T. Uehara, M. Mikuni “Sex and age dependencies of cerebral blood volume changes during cognitive activation: a multichannel near-infrared spectroscopy study”, NeuroImage 22:1715–1721 (2004)
