Open Access
Add to BibSonomyAbstract
The goal of this study was to investigate the ability of independent component analysis in the time-spectral domain to isolate physiological sources of functional near infrared spectroscopy signals. We apply independent component analysis to the broadband fNIRS data acquired on the human forehead at 650 different wavelengths between 700 nm and 950 nm. To induce cerebral oxygenation changes we use the breath holding paradigm. We found one major independent component during baseline and two major components during exercise. Each independent component corresponds to one oxy-hemoglobin and one deoxy-hemoglobin time courses. The corresponding characteristic spectra of changes in optical absorption suggested that one component represented vasodilation of cerebral arterioles while the delayed component represented the washout of deoxyhemoglobin either in cerebral capillaries and venules or in extra cerebral tissue. We found that both broadband and isolated wavelength data can produce similar independent components.
©2011 Optical Society of America
1. Introduction
Functional near infrared spectroscopy (fNIRS) aims at non-invasively measuring the effects of the brain tissue on near-infrared light propagating through the intact human head [1,2]. Various fNIRS systems use continuous-wave (CW), frequency-domain and time-resolved sources of light [1]. Although the time-resolved and frequency-domain systems have capabilities of separating absorption and scattering properties of tissues, they are more expensive and have lower signal-to-noise ratio. Therefore, CW fNIRS is considered to be most practical. CW systems use either broadband [3,4] or selected wavelength light sources [5]. While both types of sources have their advantages, a number of recent studies demonstrated benefits of broadband systems [3,4]. A drawback of all types of fNIRS systems is that the light integrates signals from all types of small blood vessels: arterioles, capillaries and venules [1,2]. Also, besides the brain the light travels through all extracerebral head tissues, whose vasomotion biases functional signals. Therefore the separation of functional signals from contamination by systemic and local physiological fluctuations is an important research topic. This problem was approached by using various signal processing methods, including blind signal separation techniques [5–8]. In particular, principal component and independent component analyses (PCA and ICA, respectively) were applied separately either to the oxy- or deoxyhemoglobin time series reconstructed from the bi-wavelength NIRS data acquired on the human head during repetitive hyper- or hypocapnia episodes [5] with the goal to filter out “surface effects”. The ICA criterion for the “surface effects” was the value of the coefficient of spatial uniformity (CSU) calculated using elements of the spatial mixing matrix [7]. In [5] the independent components with higher CSU were excluded and remaining components were recombined into the filtered oxy- and deoxyhemoglobin traces. ICA was also used in NIRS for the motion artifact removal [7–9], and to study the functional cerebral connectivity of the resting state [10].
Previously fNIRS ICA was performed only in the time-spatial domain on data acquired either at multiple source-detector separations or at different locations [6–10]. The goal of our present study was to investigate fNIRS independent components not in the time-spatial but in the time-spectral domain. In particular we wanted to understand relationships between the spectral independent components and the changes in the oxy- and deoxyhemoglobin concentrations. We show that the rows of the time-spectral ICA mixing matrix allow for a physiological interpretation in terms of the oxy- and deoxyhemoglobin concentrations. We apply a well-established Matlab version of ICA FastICA [11] to the broadband NIRS data acquired at one location on the human forehead at a single source-detector separation but at 650 different wavelengths between 700 nm and 950 nm. To induce cerebral oxygenation changes the breath holding (BH) paradigm was used because it generated blood carbon dioxide (CO2) [12]. Vasodilation induced by CO2 caused increased cerebral blood flow which resulted in changes in cerebral blood volume and oxygenation [13]. Although the individual wavelength signals were quite diverse and noisy, we found very few independent components whose origin can be interpreted based on their characteristic spectral signatures and time courses.
2. Materials and methods
2.1. Participants and experimental paradigm
Optical data was obtained from 7 healthy adult participants (1 male, 6 females), all of whom were right-handed. All participants gave informed consent before participation and the experiment has been performed according to Ryerson University 2008-003-01 Research Ethics protocol. Participants were audibly cued to perform a breath hold at the end of expiration with voluntary resumption of breathing after holding for about 15-20 seconds. The exercise was repeated six times with 40-second rest intervals. Breath holding was used in many studies (for example [3,14]) to induce hypercapnia, which typically resulted in increased concentrations of oxyhemoglobin and decreased concentrations of deoxyhemoglobin . All data set including the baseline were acquired during 6 minutes.
2.2. Data acquisition
Near-infrared light was generated by a stabilized fan-cooled AvaLight-HAL Tungsten Halogen Light Source (Avantes Inc., Broomfield, CO) with an adjustable focusing connector to maximize light coupling with the source fiber. The source fiber bundle was made of 30 Thorlabs broadband silica 400 core diameter fibers. On the probe the source fibers were
arranged circularly around the location of the detector bundle at a radius of 25 mm (Fig. 1 ). Such geometry allowed for reducing the effects of the scalp’s large blood vessels occasionally located near a light source. Light was collected using a 3-mm diameter fiber optic bundle (Sunoptic Technologies, FL). The detector bundle was connected to a QEB65000 cooled spectrometer (Ocean Optics, Dunedin, FL), which had 1044 wavelength channels in a spectral range between 650 and 1100 nm. The spectrometer output was digitized using the Spectral Suite software (Ocean Optics, Dunedin, FL). The optical probe was positioned on the left side of the forehead near the hairline to avoid sinuses. Spectra were acquired at the sampling rate of one spectrum per second resulting in 360 time points per data set. This sampling rate was selected to ensure that the instrumental noise did not affect physiological data.
2.3. Data processing: independent component analysis and general linear modeling
For 630 signals at different near-infrared wavelengths between 700 nm and 950 nm the absorbance A(λ,t) was calculated as
where Dark was the dark signal, reference was the signal measured directly through air (without tissue), and λ was the optical wavelength. Each instantaneous spectrum was smoothed using the Matlab smooth function at 15 points. No temporal smoothing was performed. The 630 x 360 matrix of 630 detrended time series of absorbance changes was fed into the FastICA algorithm [8], which used fpica Matlab function to separate the input into the N independent components (ICs):
where was the 360 x N matrix of statistically independent temporal components (ICs) [3], was the N x 630 mixing matrix, and was the matrix of residuals. The columns of were the time courses of ICs and rows of were the time-independent amplitude spectra of the corresponding components. FastICA automatically normalized the time course and the spectrum of each component such that the standard deviation (STD) of the time course was equal to one, so it had arbitrary normalized units, and the corresponding spectrum had units of the optical density. By choosing different numbers of ICs we found that two components always covered more than 90% of non-zero eigenvalues of the covariance matrix, and at the same conditions the norm of was much smaller than the norm of at each time moment. If the number of ICs N was set greater than two, the output of FastICA showed one or two non-noisy ICs (such as red traces in Fig. 1 (a,c) whose spectra spanned much greater ranges of values than the spectra of other ICs. ICs whose time courses appeared noisy in Figs. 1(a) and 1(c) also had noisy spectra (blue, green and black curves in Fig. 1(b) and blue and green curves in Fig. 1(d)). Therefore we focused on the analysis of only two components spanning the greatest ranges of values. In this approximation and for each moment of time Eq. (2) reduced to
where and were the instantaneous values of the temporal independent components, and and were the corresponding time-independent amplitude spectra. One should note that according to Eq. (3) the signs of and were arbitrary up to their simultaneous conjoined inversion. In order to homogenize the cross-subject analysis we fixated the signs of and as they appear in Fig. 1(d). This also fixated the signs of and . We assumed that, apart from the noise, the dynamic of was only due to changes in the concentrations of oxy-hemoglobin and deoxy-hemoglobin , and that it follows the modified Lambert-Beers law [1]:
where was the differential pathlength factor as a function of wavelength and r was the source detector distance. was calculated as [15] using the baseline absorbance. The baseline absorption coefficient was calculated as
where was the concentration and was the extinction coefficient of the chromophore X. The values of chromophore concentrations were taken from [16] and [17]. The concentrations of , , and water were assumed 55(micromolars), 20 , and 85%, respectively in accordance with values known from [18] for water and from [1,19] for hemoglobin.
In order to analyze the contributions of and changes to each IC, and in accordance with Eqs. (3) and (4), we performed the generalized linear model (GLM) fit (glmfit function of MATLAB) of and:
where and were the regressors of GLM, and were the parameters of the fit, and the upper indexes indicated particular ICs. It followed from Eqs. (3), (4), and (6) that for each of two ICs the changes in oxy- and deoxyhemoglobin concentrations were
respectively. Note that according to Eqs. (7) for the given k both and were synchronous in time with their time course given by . Since STDs of were equal to one, the coefficients and for each IC provided estimates of the mean magnitudes of changes in and HHb, respectively. The full and HHb changes were obtained by adding partial changes corresponding to both components:
3. Results
Figure 2 showed examples of the time courses (2a and 2c) and the corresponding spectra (2b and 2d) during rest and exercise. The main differences between rest and exercise were the number of physiological ICs (one during rest and two during exercise) and the time course of ICs (chaotic during rest and showing responses to BH during exercise). Since the shapes of the non-noisy spectra shown in Figs. 2(b,d) were qualitatively similar in all subjects we referred to the corresponding independent components as and as . The only physiological component showing non-periodic baseline hemodynamic fluctuations during rest was (identified based on the shape of its spectrum). The shape of (the spectrum of ) was topologically similar for all subjects but in some cases quantitatively changed considerably between rest and exercise. appeared only during exercises. The topological differences between and were shown in Fig. 2(d) (red and black curves, respectively).
Fig. 2 Temporal (left) and spectral (right) components during rest (a and b) and breath hold (c and d) conditions. Respective spectral and temporal components are shown by the same colors. Encircled is the IC2 spectral feature near the peak wavelength of the deoxy-hemoglobin absorption.
More examples of both component spectra were shown in Fig. 3 . Different subjects were indicated by different colors. One can see that the main differences between and were the shape of the spectrum near 760 nm (a much more pronounced spectral feature in than in near the absorption peak wavelength) and the range of values on the y axis. While the values of were all positive, the values of were either all negative or negative at shorter wavelengths and positive at longer wavelengths. We note that although the signs of were arbitrary up to their conjoined inversion with , one can see from Fig. 3 that the inversions of signs could not convert spectra of one type into another type. Therefore, we used the topological differences between the shapes of and for the initial identifications of and for each subject. The physiology-related quantitative differences between and were revealed by the results of GLM fitting shown in Table 1 . From Table 1 one can see that the main difference between two components was in the relative magnitude of the oxy-hemoglobin changes, represented by compared to that of the deoxy-hemoglobin, represented by .
Table 1. Results of the fit of using GLM in Eq. (5)
Since was much greater than, the contribution of the deoxy-hemoglobin changes to was much smaller than the contribution of the oxy-hemoglobin changes. For on average both oxy- and deoxy-hemoglobin changes were approximately equal by magnitude and, on average, opposite by sign. Figure 4 shows examples of the time courses of and for three subjects during the repetitive BH exercises. In all panels of Fig. 3 the increases in corresponded to increases in . The direction of change in depended on the sign of in Table 1. The increases in in Figs. 3 (b-d) corresponded to the increases in concentration and to the decreases in . One can see that for different subjects the responses to BHs (Fig. 4 (a)) were all quite synchronous. The responses were always delayed compared to the responses (Figs. 4 (b, c, d). Delays between the peaks of and were different for different subjects, varying from 4s to 15s. We have verified that the correlation coefficient between and was always zero.
Fig. 4 (a) Time courses of IC1 for three different subjects during BH. (b-d) time courses of IC1 (red curve) and IC2 (blue curve) for three sample subjects. The components were extracted from the broadband data sets. Vertical bars showed BH periods. Note that the time courses were normalized to have their STDs equal to one.
Figure 5 shows the full and traces computed according to Eqs. (8), and the total-hemoglobin changes for the same data set as Fig. 4b. The scaled and were also shown. The figure demonstrates that and were well correlated with , while was better correlated with . However, these features did not necessarily mean that represented mainly in the brain and represented the cerebral . If the actual sources of and were located in different compartments of the vascular bed, the similarities of to and of to in Fig. 6 could result from the purely numeric fact that the values of were much greater than all other coefficients and conjointly the values of were much smaller than all other coefficients in Table 1.
Fig. 5 The full ΔHbO2 (t) and ΔHHb(t) traces computed according to Eqs. (8), and the total-hemoglobin changes. T1(t) and T2(t) are also shown by dashed and dotted lines, respectively.
Fig. 6 Comparison of temporal components IC1 (a) and IC2 (b), obtained from broadband data (solid line) and using two narrow wavebands n 760 and 900 nm (dashed line).
In order to investigate whether and could be obtained using a small number of wavelengths, we created five signals by averaging ten adjacent wavelength signals around each of 760 nm, 800 nm, 830 nm, 860 nm, and 900 nm wavelengths. Then and were obtained by applying FastICA to the pairs of signals, one at 760 nm and another either at 800 nm, or at 830 nm, or at 860 nm, or 900 nm. The comparisons of the time courses of ICs obtained using the full spectrum and a pair of wavelengths (760 nm and 900 nm) were presented in Fig. 5. One can see that the time courses obtained using the full and reduced spectra were similar and the main difference was a much noisier obtained using two wavelengths. Figure 6 compares ICs obtained using the full spectrum and pairs of wavelengths in terms of STDs of the differences between the corresponding temporal ICs. We note that the values of standard deviations in Fig. 6 should be also compared with the standard deviations of and , which were both equal to one due to the normalization. From Fig. 6 one can conclude that (a) for all pairs of wavelengths the differences in from the full spectrum were much larger than for , (b) for the difference was minimal at 860 nm, and (c) for the difference was minimal at 830 nm. Also from Fig. 6 one can see that the wavelength of 830 nm allowed for about 15-20% improvement compared to 800 nm and 900 nm. The effect of the wavelength tuning on was very modest. Also we found that using combined signals at all five wavelengths (760 nm, 800 nm, 830 nm, 860 nm, and 900 nm) gave almost the same accuracy of as the pair of 760 nm and 860 nm, and the same accuracy of as the pair of 760 nm and 830 nm (Fig. 7 ).
Fig. 7 Standard deviations of IC time courses computed using the full spectrum relative to those obtained using two wavelengths: one at 760 nm and another at a wavelength shown.
4. Discussion
Unlike other groups who applied ICA to fNIRS data in the time-spatial domain, we performed ICA on the broadband fNIRS data in the time-spectral domain. While others interpreted the spatial ICs based on their spatial uniformity [6,7] the spectral approach allows for the component interpretation based on the and dynamic contributions to each component. Generally, ICA separates the n - channel input signals into N statistically independent components. The statistical independence means that the joint probability of components to have certain values equals to the product of the individual probabilities. In particular, for time variables having zero means the statistical independence also results in their temporal uncorrelatedness. Although multichannel signals of any nature can be separated into ICs, they should not necessarily be associated with real physiological processes of different physical origin. However, the fact discovered in this study that two major temporal ICs had their optical spectra and relative time courses similar for many subjects could indicate that the physiological sources of these components were different. Our results allowed for two different interpretations. First was that and corresponded to the cerebral oxy- and deoxy-hemoglobin changes, respectively (as suggested by Fig. 5). Another possible interpretation was that both and corresponded to both oxy- and deoxy-hemoglobin changes but in different parts of the vascular bed. From this perspective, represented synchronous changes in both hemoglobin species in opposite directions, while the small coefficient for showed that represented mostly changes in the volume of highly oxygenated blood. The facts that existed during the baseline and led the responses to BHs and that had much more contribution from than from suggested that could be related to the spontaneous or CO2-induced dilation of arterioles in the brain. The delayed and opposite - responses represented by could be related either to the washout of in capillaries and small veins, or to the exercise-related hemodynamics in extracerebral tissue. Note that the increases in in Figs. 3 (b-d) corresponded to the increases in concentration and to the decreases in . The spectral features of and and their related physiological meaning would be very difficult to recognize using a few selected wavelengths instead of our broadband system, which allowed us to identify spectral differences shown in Fig. 3. To some extent these issues can be investigated using fMRI. Blood oxygen level dependent (BOLD) fMRI data obtained by different groups during hypercapnia show BOLD responses rather similar to the traces in Fig. 4 (a) than to time course in Figs. 4 (b - d) [3,20]. This supports the hypothesis that was related to the cerebral hemodynamics while resulted from a different source. However, one should note that the exact source of fMRI signal depends on the strength of magnetic field and other factors [21]. Another possible approach to further study the origin of discovered ICs, in particular to test if one or both components originated from the scalp is to use a combination of multi-wavelength measurement channels with small and large source-detector separations.
We found that the responses from 630 spectral channels to BHs can be represented by only two ICs. Such a dramatic reduction of the data dimension creates a question of whether instead of the entire spectrum one can use a few isolated wavelengths. Since the signal-to-noise ratio is an important issue of all types of measurements, a complete answer to this question requires an additional study comparing data obtained using the broadband and multi-wavelength devices. However, the results presented in Fig. 6 showed that the time courses of ICs obtained using the full broadband spectrum and two wavelengths at 760 nm and 900 nm were similar. We stress that ICA applied to the two-wavelength data returned two ICs each corresponding to both and changes, and a 2 x 2 mixing matrix, whose four elements together with the extinction coefficients and can be used to calculate contributions of both and to each component. Therefore, if different ICs indeed have different and typical physiological sources, ICA can be used to separate contributions of those sources to data acquired at one location and single source-detector separation even if the number of wavelengths is small. This approach can be combined with the multi-distance approach (such as in [6]) to better separate extracerebral and cerebral signals. The results of our investigation of the reproducibility of the time courses of the broadband ICs using few selected wavelengths showed that was much less affected by the spectrum reduction than and that the tuning of isolated wavelengths or increasing of their number caused only modest improvements.
5. Conclusion
By applying independent component analysis to broadband fNIRS data we found one major IC during baseline and two major temporal ICs, one leading and another delayed, during repetitive breath holds. The corresponding characteristic spectra of changes in optical absorption suggested that the leading component represented vasodilation of cerebral arterioles while the delayed component represented the washout of deoxyhemoglobin in either cerebral venules or in extracerebral tissue.
References and links
1. D. A. Benaron, S. R. Hintz, A. Villringer, D. Boas, A. Kleinschmidt, J. Frahm, C. Hirth, H. Obrig, J. C. van Houten, E. L. Kermit, W. F. Cheong, and D. K. Stevenson, “Noninvasive functional imaging of human brain using light,” J. Cereb. Blood Flow Metab. 20(3), 469–477 (2000). [CrossRef] [PubMed]
2. F. Irani, S. M. Platek, S. Bunce, A. C. Ruocco, and D. Chute, “Functional near infrared spectroscopy (fNIRS): an emerging neuroimaging technology with important applications for the study of brain disorders,” Clin. Neuropsychol. 21(1), 9–37 (2007). [CrossRef] [PubMed]
3. I. Schelkanova and V. Toronov, “Optimal quantitation of the cerebral hemodynamic response in functional near-infrared spectroscopy,” Opt. Express 18(18), 19386–19395 (2010). [CrossRef] [PubMed]
4. H. Dehghani, F. Leblond, B. W. Pogue, and F. Chauchard, “Application of spectral derivative data in visible and near-infrared spectroscopy,” Phys. Med. Biol. 55(12), 3381–3399 (2010). [CrossRef] [PubMed]
5. J. Virtanen, T. Noponen, and P. Meriläinen, “Comparison of principal and independent component analysis in removing extracerebral interference from near-infrared spectroscopy signals,” J. Biomed. Opt. 14(5), 054032 (2009). [CrossRef] [PubMed]
6. A. Hyvarinen, J. Karhunen, and E. Oja, Independent Component Analysis (Wiley, New York, 2001).
7. S. Kohno, I. Miyai, A. Seiyama, I. Oda, A. Ishikawa, S. Tsuneishi, T. Amita, and K. Shimizu, “Removal of the skin blood flow artifact in functional near-infrared spectroscopic imaging data through independent component analysis,” J. Biomed. Opt. 12(6), 062111 (2007). [CrossRef] [PubMed]
8. F. C. Robertson, T. S. Douglas, and E. M. Meintjes, “Motion artifact removal for functional near infrared spectroscopy: a comparison of methods,” IEEE Trans. Biomed. Eng. 57(6), 1377–1387 (2010). [CrossRef] [PubMed]
9. A. V. Medvedev, J. M. Kainerstorfer, S. V. Borisov, A. H. Gandjbakhche, and J. Vanmeter, ““Seeing” electroencephalogram through the skull: imaging prefrontal cortex with fast optical signal,” J. Biomed. Opt. 15(6), 061702 (2010). [CrossRef] [PubMed]
10. H. Zhang, Y. J. Zhang, C. M. Lu, S. Y. Ma, Y. F. Zang, and C. Z. Zhu, “Functional connectivity as revealed by independent component analysis of resting-state fNIRS measurements,” Neuroimage 51(3), 1150–1161 (2010). [CrossRef] [PubMed]
11. C. Hesse and C. James, “The fastica algorithm with spatial constraints,” IEEE Signal Process. Lett. 12(11), 792–795 (2005). [CrossRef]
12. T. Len, J. Neary, G. Asmundson, D. Goodman, B. Bjornson, and Y. Bhambhani, “Cerebrovascular reactivity impairment after sport-induced concussion,” Med. Sci. Sports Exerc. 43(12), 2241–2248 (2011). [CrossRef]
13. R. L. Grubb Jr, M. E. Raichle, J. O. Eichling, and M. M. Ter-Pogossian, “The effects of changes in PaCO2 on cerebral blood volume, blood flow, and vascular mean transit time,” Stroke 5(5), 630–639 (1974). [CrossRef] [PubMed]
14. U. E. Emir, C. Ozturk, and A. Akin, “Multimodal investigation of fMRI and fNIRS derived breath hold BOLD signals with an expanded balloon model,” Physiol. Meas. 29(1), 49–63 (2008). [CrossRef] [PubMed]
15. A. Sassaroli and S. Fantini, “Comment on the modified Beer-Lambert law for scattering media,” Phys. Med. Biol. 49(14), N255–N257 (2004). [CrossRef] [PubMed]
16. W. B. Gratzer, Medical Research Council Laboratories, Holly Hill, London, personal communication.
17. N. Kollias, Wellman Laboratories, Harvard Medical School, Boston, MA, personal communication.
18. H. Q. Woodard and D. R. White, “The composition of body tissues,” Br. J. Radiol. 59(708), 1209–1218 (1986). [CrossRef] [PubMed]
19. J. H. Choi, M. Wolf, V. Toronov, U. Wolf, C. Polzonetti, D. Hueber, L. P. Safonova, R. Gupta, A. Michalos, W. Mantulin, and E. Gratton, “Noninvasive determination of the optical properties of adult brain: near-infrared spectroscopy approach,” J. Biomed. Opt. 9(1), 221–229 (2004). [CrossRef] [PubMed]
20. B. J. MacIntosh, L. M. Klassen, and R. S. Menon, “Transient hemodynamics during a breath hold challenge in a two part functional imaging study with simultaneous near-infrared spectroscopy in adult humans,” Neuroimage 20(2), 1246–1252 (2003). [CrossRef] [PubMed]
21. B. A. Poser, E. van Mierlo, and D. G. Norris, “Exploring the post-stimulus undershoot with spin-echo fMRI: implications for models of neurovascular response,” Hum. Brain Mapp. 32(1), 141–153 (2011). [CrossRef] [PubMed]
