Multiregional functional near-infrared spectroscopy reveals globally symmetrical and frequency-specific patterns of superficial interference

Linear regression with short source-detector separation channels (S-channels) as references is an efficient way to overcome significant physiological interference from the superficial layer for functional nearinfrared spectroscopy (fNIRS). However, the co-located configuration of Schannels and long source-detector separation channels (L-channels) is difficult to achieve in practice. In this study, we recorded superficial interference with S-channels in multiple scalp regions. We found that superficial interference has overall frequency-specific and globally symmetrical patterns. The performance of linear regression is also dependent on these patterns, indicating the possibility of simplifying the Schannel configurations for multiregional fNIRS imaging. ©2015 Optical Society of America OCIS codes: (170.2655) Functional monitoring and imaging; (170.3880) Medical and biological imaging; (170.5380) Physiology; (170.6960) Tomography References and links 1. S. C. Bunce, M. Izzetoglu, K. Izzetoglu, B. Onaral, and K. Pourrezaei, “Functional near-infrared spectroscopy,” IEEE Eng. Med. Biol. Mag. 25(4), 54–62 (2006). 2. Y. Hoshi, “Functional near-infrared spectroscopy: current status and future prospects,” J. Biomed. Opt. 12(6), 062106 (2007). 3. 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). 4. Y. Minagawa-Kawai, H. van der Lely, F. Ramus, Y. Sato, R. Mazuka, and E. Dupoux, “Optical brain imaging reveals general auditory and language-specific processing in early infant development,” Cereb. Cortex 21(2), 254–261 (2011). 5. R. J. Cooper, J. C. Hebden, H. O’Reilly, S. Mitra, A. W. Michell, N. L. Everdell, A. P. Gibson, and T. Austin, “Transient haemodynamic events in neurologically compromised infants: a simultaneous EEG and diffuse optical imaging study,” Neuroimage 55(4), 1610–1616 (2011). 6. J. Huang, F. Wang, Y. Ding, H. Niu, F. Tian, H. Liu, and Y. Song, “Predicting N2pc from anticipatory HbO activity during sustained visuospatial attention: A concurrent fNIRS-ERP study,” Neuroimage 113, 225–234 (2015). 7. T. J. Huppert, R. D. Hoge, S. G. Diamond, M. A. Franceschini, and D. A. Boas, “A temporal comparison of BOLD, ASL, and NIRS hemodynamic responses to motor stimuli in adult humans,” Neuroimage 29(2), 368–382 (2006). 8. R. D. Hoge, M. A. Franceschini, R. J. Covolan, T. Huppert, J. B. Mandeville, and D. A. Boas, “Simultaneous recording of task-induced changes in blood oxygenation, volume, and flow using diffuse optical imaging and arterial spin-labeling MRI,” Neuroimage 25(3), 701–707 (2005). 9. D. A. Boas, A. M. Dale, and M. A. Franceschini, “Diffuse optical imaging of brain activation: Approaches to optimizing image sensitivity, resolution, and accuracy,” Neuroimage 23(Suppl 1), S275–S288 (2004). 10. R. C. Mesquita, M. A. Franceschini, and D. A. Boas, “Resting state functional connectivity of the whole head with near-infrared spectroscopy,” Biomed. Opt. Express 1(1), 324–336 (2010). 11. C. M. Lu, Y. J. Zhang, B. B. Biswal, Y. F. Zang, D. L. Peng, and C. Z. Zhu, “Use of fNIRS to assess resting #236645 Received 23 Mar 2015; revised 25 May 2015; accepted 16 Jun 2015; published 8 Jul 2015 (C) 2015 OSA 1 Aug 2015 | Vol. 6, No. 8 | DOI:10.1364/BOE.6.002786 | BIOMEDICAL OPTICS EXPRESS 2786 state functional connectivity,” J. Neurosci. Methods 186(2), 242–249 (2010). 12. B. R. White, A. Z. Snyder, A. L. Cohen, S. E. Petersen, M. E. Raichle, B. L. Schlaggar, and J. P. Culver, “Resting-state functional connectivity in the human brain revealed with diffuse optical tomography,” Neuroimage 47(1), 148–156 (2009). 13. Y. J. Zhang, C. M. Lu, B. B. Biswal, Y. F. Zang, D. L. Peng, and C. Z. Zhu, “Detecting resting-state functional connectivity in the language system using functional near-infrared spectroscopy,” J. Biomed. Opt. 15(4), 047003 (2010). 14. Q. Zhang, E. N. Brown, and G. E. Strangman, “Adaptive filtering to reduce global interference in evoked brain activity detection: a human subject case study,” J. Biomed. Opt. 12(6), 064009 (2007). 15. F. Tian, H. Niu, B. Khan, G. Alexandrakis, K. Behbehani, and H. Liu, “Enhanced Functional Brain Imaging by Using Adaptive Filtering and a Depth Compensation Algorithm in Diffuse Optical Tomography,” IEEE Trans. Med. Imaging 30(6), 1239–1251 (2011). 16. Q. Zhang, G. E. Strangman, and G. Ganis, “Adaptive filtering to reduce global interference in non-invasive NIRS measures of brain activation: how well and when does it work?” Neuroimage 45(3), 788–794 (2009). 17. R. B. Saager and A. J. Berger, “Direct characterization and removal of interfering absorption trends in two-layer turbid media,” J. Opt. Soc. Am. A 22(9), 1874–1882 (2005). 18. B. W. Zeff, B. R. White, H. Dehghani, B. L. Schlaggar, and J. P. Culver, “Retinotopic mapping of adult human visual cortex with high-density diffuse optical tomography,” Proc. Natl. Acad. Sci. U.S.A. 104(29), 12169– 12174 (2007). 19. N. M. Gregg, B. R. White, B. W. Zeff, A. J. Berger, and J. P. Culver, “Brain specificity of diffuse optical imaging: improvements from superficial signal regression and tomography,” Front Neuroenergetics 2, 14 (2010). 20. R. Saager and A. Berger, “Measurement of layer-like hemodynamic trends in scalp and cortex: implications for physiological baseline suppression in functional near-infrared spectroscopy,” J. Biomed. Opt. 13(3), 034017 (2008). 21. S. Umeyama and T. Yamada, “Monte Carlo study of global interference cancellation by multidistance measurement of near-infrared spectroscopy,” J. Biomed. Opt. 14(6), 064025 (2009). 22. T. Yamada, S. Umeyama, and K. Matsuda, “Multidistance probe arrangement to eliminate artifacts in functional near-infrared spectroscopy,” J. Biomed. Opt. 14(6), 064034 (2009). 23. L. Gagnon, K. Perdue, D. N. Greve, D. Goldenholz, G. Kaskhedikar, and D. A. Boas, “Improved recovery of the hemodynamic response in diffuse optical imaging using short optode separations and state-space modeling,” Neuroimage 56(3), 1362–1371 (2011). 24. M. A. Franceschini, D. K. Joseph, T. J. Huppert, S. G. Diamond, and D. A. Boas, “Diffuse optical imaging of the whole head,” J. Biomed. Opt. 11(5), 054007 (2006). 25. Y. Zhang, D. H. Brooks, M. A. Franceschini, and D. A. Boas, “Eigenvector-based spatial filtering for reduction of physiological interference in diffuse optical imaging,” J. Biomed. Opt. 10(1), 011014 (2005). 26. 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). 27. L. Gagnon, R. J. Cooper, M. A. Yücel, K. L. Perdue, D. N. Greve, and D. A. Boas, “Short separation channel location impacts the performance of short channel regression in NIRS,” Neuroimage 59(3), 2518–2528 (2012). 28. L. Gagnon, M. A. Yücel, D. A. Boas, and R. J. Cooper, “Further improvement in reducing superficial contamination in NIRS using double short separation measurements,” Neuroimage 85(Pt 1), 127–135 (2014). 29. E. Kirilina, A. Jelzow, A. Heine, M. Niessing, H. Wabnitz, R. Brühl, B. Ittermann, A. M. Jacobs, and I. Tachtsidis, “The physiological origin of task-evoked systemic artefacts in functional near infrared spectroscopy,” Neuroimage 61(1), 70–81 (2012). 30. S. B. Erdoğan, M. A. Yücel, and A. Akın, “Analysis of task-evoked systemic interference in fNIRS measurements: insights from fMRI,” Neuroimage 87, 490–504 (2014). 31. A. T. Eggebrecht, S. L. Ferradal, A. Robichaux-Viehoever, M. S. Hassanpour, H. Dehghani, A. Z. Snyder, T. Hershey, and J. P. Culver, “Mapping distributed brain function and networks with diffuse optical tomography,” Nat. Photonics 8(6), 448–454 (2014). 32. G. H. Klem, H. O. Lüders, H. H. Jasper, and C. Elger; The International Federation of Clinical Neurophysiology, “The ten-twenty electrode system of the International Federation,” Electroencephalogr. Clin. Neurophysiol. Suppl. 52, 3–6 (1999). 33. M. Cope and D. T. Delpy, “System for long-term measurement of cerebral blood and tissue oxygenation on newborn infants by near infra-red transillumination,” Med. Biol. Eng. Comput. 26(3), 289–294 (1988). 34. B. R. White, S. M. Liao, S. L. Ferradal, T. E. Inder, and J. P. Culver, “Bedside optical imaging of occipital resting-state functional connectivity in neonates,” Neuroimage 59(3), 2529–2538 (2012). 35. P. D. Welch, “Use of fast Fourier transform for estimation of power spectra a method based on time averaging over short modified periodograms,” IEEE Trans. Acoust. Speech 15, 70–73 (1967). 36. D. P. Auer, “Spontaneous low-frequency blood oxygenation level-dependent fluctuations and functional connectivity analysis of the ‘resting’ brain,” Magn. Reson. Imaging 26(7), 1055–1064 (2008). 37. B. Biswal, F. Z. Yetkin, V. M. Haughton, and J. S. Hyde, “Functional connectivity in the motor cortex of resting human brain using echo-planar MRI,” Magn. Reson. Med. 34(4), 537–541 (1995). 38. M. D. Fox and M. E. Raichle, “Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging,” Nat. Rev. Neurosci. 8(9), 700–711 (2007). 39. F. Tian and H. Liu, “Depth-compensated diffuse optical tomography enhanced by general linear model analysis and an anatomical atlas of human head,” Neuroimage 85(Pt 1), 166–180 (2014). #236645 Received 23 Mar 2015; revised 25 May 2015; accepted 16 Jun 2015; published 8 Jul 2015 (C) 2015 OSA 1 Aug 2015 | Vol. 6, No. 8 | DOI:10.1364/BOE.6.002786 | BIOMEDICAL OPTICS EXPRESS 2787 40. Y. Tong, L. M. Hocke, L. D. Nickerson, S. C. Licata, K. P. Lindsey, and B. Frederick, “Evaluating the effects of systemic low frequency oscillations measured in the periphery on the independent component analysis results of resting state networks,” Neuroimage 76, 202–215 (2013).


Introduction
Functional near-infrared spectroscopy (fNIRS) is an emerging and promising brain imaging technique.It has become an important complement to functional Magnetic Resonance Imaging (fMRI) for measuring hemodynamic brain activities with several unique advantages: First, fNIRS is low-cost, salient and portable, making it possible for long-term monitoring and repeated measurements of brain activities in various situations [1][2][3].Such properties also make fNIRS suitable for the pediatric population and patients who cannot be scanned in fMRI [4,5].Second, the compatibility of fNIRS with other electrical or magnetic systems, such as electroencephalography (EEG) [5,6], favors multi-modal imaging of brain activities that is important for investigating the brain at multiple spatial and temporal scales simultaneously.Finally, fNIRS measures changes in oxygenated hemoglobin (HbO) and deoxygenated hemoglobin (HbR), delivering additional information concerning the cerebral metabolic state compared to the blood-oxygen-level dependent (BOLD) signal in fMRI [7,8].
In fNIRS imaging, the near-infrared light is emitted by a source that is placed on the scalp.The light travels through the superficial layer (that is, the scalp and skull), then reaches the cortical tissues, and finally returns to a detector that is on the same side of the scalp.Therefore, the measured signal is sensitive not only to the hemodynamic activities in the brain, but also to the vascular fluctuations in the superficial layer.The latter mainly refers to a number of systemic physiological fluctuations arising from the cardiac pulsations (around 1 to 2 Hz), respiration (around 0.2 to 0.4 Hz), Mayer waves (around 0.1 Hz), and other very lowfrequency fluctuations (0.01 to 0.05 Hz) [9].These physiological fluctuations may induce significant interference to the desired signal from the cerebral tissues in both the resting state [10][11][12][13] and the task-evoked state [14][15][16].Thus, the reduction of superficial interference is one of the main challenges in fNIRS.
Linear regression of superficial interference is an effective method to address this issue [17].In addition to long source-detector separation channels (L-channels where sourcedetector separation is about 3 to 4 cm) that are sensitive to both the superficial layer and the brain, some short source-detector separation channels (S-channels where source-detector separation is about 0.5 to 1.5 cm) can be used to measure the physiological fluctuations within the superficial layer only [12,18,19].By assuming that superficial interference measured by the S-channels is homogeneous to the corresponding components in the Lchannels [20], the components of superficial interference in the L-channels can be regressed out by using the S-channel as a reference.
Although the effectiveness of this S-channel-based regression method has been proved empirically [12,[21][22][23], how to arrange the S-channels along with the L-channels is still debatable.Early studies conducting S-channel-based regression relied on the assumption that superficial interference is globally homogenous across the scalp.Single or a few S-channels were used for all L-channels [17,21,22].This assumption was justified by some indirect evidence.For example, by correlating the fNIRS signals from several L-channels with the peripheral physiological signals measured by auxiliary devices, Franceschini et al. found that the L-channel fNIRS signals over the head were coherent with the peripheral physiological fluctuations [24].A few studies using principle component analysis (PCA) [25] or independent component analysis (ICA) [26] also found that the physiological components derived from L-channel measurements were globally distributed.However, these studies were not focused on the physiological fluctuations in the superficial layer.In fact, the physiological components in the L-channel signals in the above studies originated from both the superficial layer and the brain.
On the other hand, a few studies have suggested that superficial interference is spatially inhomogeneous across the scalp.Gagnon and his associates [27] reported that, the correlation between S-and L-channels and the denoising performance of S-channel-based regression decreased as the relative distance between the S-channels and L-channels increased in the left motor cortex.They suggested that the relative distance between L-channel and its corresponding S-channel should be less than 1.5 cm.In a later study, they further suggested that ideally two S-channels should be used for one L-channel, one at the source optode and the other at the detector optode of the L-channel [28].A study by Kirilina et al. [29] showed that task-evoked systemic artifacts in the forehead were not homogeneously distributed, but were rather localized in the scalp draining veins.Furthermore, based on Monte Carlo simulations, Erdogan et al. reported that the de-noising performance of linear regression was greatest when using both the global signals and the local S-channel signals as regressors for L-channel signals in the forehead, implying that superficial interference is spatially inhomogeneous at that location.Therefore, they concluded that maximizing the overlap between the S-and L-channels is of great importance [30].However, these pioneering studies were conducted within a limited area (the left motor cortex or the forehead).To the best of our knowledge, it is still unknown whether the spatial inhomogeneity of superficial interference affects the entire human head In recent years, multiregional fNIRS imaging has become increasingly popular to map distributed brain function and networks [31].However, the application of S-channel-based regression in multiregional fNIRS imaging has been hindered due to currently limited knowledge on the spatial similarity of superficial interference.The aim of the present study was to explore the spatial similarity of superficial interference at multiple scalp regions over the head using data acquired in a resting state and to get preliminary suggestions for configuring S-channels in multiregional fNIRS studies.Given that superficial interference has several physiological origins, that feature was also characterized in four frequency bands that are related to the cardiac pulsations, respiration, Mayer waves and very low-frequency fluctuations.Then, S-channel-based regression on data acquired during a motor task was conducted to prove the usefulness of the identified spatial similarity of superficial interference in configuring the S-channels for multiregional fNIRS imaging.

Subjects
Thirteen healthy adults were recruited for two separate experiments.Informed consents were obtained before the experiments according to the procedure approved by the Institutional Review Board at the State Key Laboratory of Cognitive Neuroscience and Learning, Beijing Normal University.

Experiments
In the first experiment, superficial interference was measured in the resting state using multiple S-channels to investigate spatial similarity across different scalp regions.The configuration of fNIRS probe is shown in Fig. 1(A).It consisted of 8 sources and 24 detectors.The resultant 24 S-channels (source-detector separation was 1.5 cm) were located in four scalp regions that were symmetrical between the two cerebral hemispheres: the anterior, middle, lateral and posterior regions.According to the international 10-20 system [32], the S-channels in the anterior region were along the Fp1-Fp2 line, the S-channels in the middle and lateral regions were along the T3-C3-C4-T4 line, and the S-channels in the posterior region were along the O1-O2 line.
Eight subjects (6 males and 2 females, age = 25.9 ± 4.1 years, all right-handed) participated in the first experiment.The subjects were instructed to sit still for 15 minutes during data collection with eyes closed and body relaxed.In the second experiment, brain activities during a resting state followed by a motor task session were measured with L-channels on the bilateral sensorimotor cortices and S-channels spread over the bilateral anterior, middle and posterior regions of the head.The fNIRS probe consisted of 10 sources and 24 detectors as shown in Fig. 2(A).According to the international 10-20 system [32], the 20 L-channels (source-detector separation was 3 cm or 3.35 cm) were positioned around C3 and C4, respectively.These L-channels mainly covered the bilateral sensorimotor cortices.The four S-channels in the middle region were placed symmetrically around C3 and C4, too.The four S-channels in the anterior region were placed along the Fp1-Fp2 line, and the four S-channels in the posterior region were placed along the O1-O2 line.Fig. 2. Geometry of fNIRS probe and types of S-channels (S-ch) in the second experiment.All the L-channels were located around the bilateral sensorimotor cortices.Using one L-channel (thick red line) as an example, all S-channels could be categorized into four types according to their locations: local S-channels (orange, in the same sensorimotor region as the selected Lchannel), symmetrical-middle S-channels (black, in the sensorimotor region symmetrical to the selected L-channel), posterior S-channels (blue), and anterior S-channels (green).previously.During the resting-state session, the subjects were instructed to sit still for about 15 minutes with eyes closed and body relaxed, just as they did in the first experiment.The subsequent motor task session was in block design: In each block, each subject tapped the index finger of both hands in a random sequence according to visual cues presented on a computer screen.Specifically, a 4-letter cue consisting of "F" and "J" (e.g., "FFJJ" or "FJJF") was presented for 2 seconds.The subjects were instructed to press the keyboard buttons "F" and "J" as quickly and correctly as possible by using their left and right index fingers, respectively.The same cue was repeated for 10-15 times (i.e., 20-30 seconds) randomly within each block, followed by a resting period that was in equal duration to the finger tapping period.The total length of the motor task session was about 6 minutes by including seven tapping-resting blocks.
In both experiments, the CW6 system from TechEn, Inc. (Medford, MA, USA) [24] was used to collect data from the subjects.The absorption of near infrared light was measured at two wavelengths of 690 nm and 830 nm with a sampling rate of 25 Hz.Custom engineering was performed in construction of the fNIRS probes.Thin polyethylene strips were used to hold the optodes in the desired source-detector separations.Specifically, three polyethylene strips were made for each experiment: one was for the measurements in the bilateral anterior regions along the Fp1-Fp2 line; one was for the measurements in the bilateral middle-lateral regions along the T3-C3-C4-T4 line, and one was for measurements in the bilateral posterior regions along the O1-O2 line.The polyethylene strips were connected with self-grip Velcro and the separations between them were carefully adjusted according to the international 10-20 measurements [32] of each individual's head.

Data preprocessing
Prior to analysis, five channels in the first experiments (among eight subjects × 24 channels/subject = 192 channels) and 30 channels in the second experiments (among nine subjects × 32 channels/subject = 288 channels) were excluded due to bad channel warnings from the CW6 system.Such warnings indicated poor or inadequate optical signals that usually resulted from hair absorption and/or imperfect contact of optodes to the scalp.
For the qualified channels, changes in HbO and HbR concentrations were calculated based on the modified Beer-Lambert Law [33].Two differential pathlength factors were 6.8 at 690 nm and 5.8 at 830 nm.Then, the channel-wise changes of HbO and HbR were detrended by removing the 1st-and 2nd-order polynomial drifts.To further reduce motion artifacts and large instrumental noise, channels with high standard deviations of HbO and/or HbR time series were removed from any further analysis.The threshold was set as the top 7.5% of the standard deviations of the time courses among all the channels and all the subjects (that is, 14 channels in the first experiment and 19 channels in the second experiment were removed) [12,34].Using Butterworth filters, the remaining data were low-pass filtered at 2 Hz to remove the noise at ultra-high frequencies and then high-pass filtered at 0.01 Hz to remove residual long-term drifts.This data set will be referred to as "whole-frequency data" hereafter.

Identifying the spatial magnitude patterns of superficial interference
Based on the whole-frequency data from the first experiment, we assessed the power of superficial interference in the four scalp regions (that is, the anterior, lateral, middle and posterior regions).First, we computed the power spectral densities of all S-channels using Welch's method with a Hamming window [35].Then, the power of the four typical physiological fluctuations was extracted as a sum in their respective frequency bands: 0.01-0.05Hz for the very low-frequency fluctuations, 0.07-0.15Hz for the Mayer waves, 0.2-0.4Hz for the respiratory waves, and 0.8-2 Hz for the cardiac pulsations [24].Finally, we calculated the mean power of all S-channels within each of the four scalp regions.

Identifying the spatial similarity patterns of superficial interference
For the first experiment, the cross-correlation function between every two paired S-channels was calculated from both the HbO and HbR time series of each subject.The time lags ranged from −50 to 50 seconds with an interval of 0.04 second.The degree of similarity between every two paired S-channels was defined as the maximum of their cross-correlation coefficients.Then, we evaluated the similarities of S-channel pairs at two different spatial levels: the local level (the two paired channels were on the same hemisphere and within one of the four regions: the anterior, middle, lateral and posterior regions) and the global level (the two paired channels were on different hemispheres and/or in different regions).At the local level, a total of 24 S-channel pairs within each of the four scalp regions (referred to as loS-ch, shown in Fig. 1(B)) were analyzed.At the global level, a total of 252 S-channel pairs across different hemispheres and/or regions were analyzed.According to their relative locations, as shown in Fig. 1(C), the global S-channel pairs were further categorized into three types: the symmetrical, ipsilateral, and contralateral-asymmetrical pairs.Specifically, there were 36 symmetrical pairs (syS-ch), 108 ipsilateral pairs (ipS-ch), and 108 contralateralasymmetrical pairs (caS-ch).Finally, we compared the similarity values among the local pairs and three types of global pairs across all the subjects using the two-tailed t test.
This entire procedure was also repeated after band-pass filtering of the whole-frequency data to extract the four typical physiological fluctuations in different frequency bands.Specifically, the whole-frequency data was band-pass filtered between 0.01 and 0.05 Hz to extract the very low-frequency fluctuations, between 0.07 and 0.15 Hz to extract the Mayer waves, between 0.2 and 0.4 Hz to extract the respiratory waves, and between 0.8 and 2 Hz to extract the cardiac pulsations [24].Then, we calculated the frequency-specific crosscorrelation coefficients between every two S-channels in the same way as we did on the whole-frequency data.Finally, we compared the similarity values of the four typical physiological fluctuations among the local pairs and three types of global pairs across all the subjects using the two-tailed t test.
For the similarities described above, values of ≥0.9 were considered very high; 0.7 to 0.9 were considered high; 0.5 to 0.7 were considered moderate, and <0.5 were considered low.

Evaluating the denoising performance of linear regression with respect to the spatial similarity patterns of superficial interference
Based on the whole-frequency data collected during the motor task in the second experiment, removal of superficial interference from the L-channels was conducted by using the Schannels at multiple scalp regions as regressors.Specifically, linear regression was performed according to the formula below: where Y represents the signal from an L-channel; X represents the signal from an S-channel with a suitable time lag [17]; and ε represents the residual of linear regression that is also referred to as the denoised result of the L-channel.In this study, the time lag between Schannel and L-channel was set as the one showing the maximum cross-correlation coefficient between the two time courses.According to the relative locations of each S-channel to an L-channel, all S-channel regressors were categorized into four types: (1) the S-channels adjacent to each L-channel (≤ 5 cm, referred to as Local S-ch regressors), as illustrated by the orange circles in Fig. 2 with the thick red line representing the interested L-channel; (2) the S-channels at the middle region that were symmetrical to the interested L-channel (referred to as Sym-Mid S-ch regressors), as illustrated by the black circles in the same figure ; (3) the S-channels at the posterior regions (referred to as Posterior S-ch regressors), as illustrated by the blue circles in the same figure; and (4) the S-channels at the anterior regions (referred to as Anterior S-ch regressors), as illustrated by the green circles in the same figure.Additionally, we also used the global mean of all S-channels over the head as a single regressor.
Because the results derived from the real experiment had a lack of ground truth in terms of task-evoked responses, we referred to the results with local regressors as a standard because of their excellent denoising performance, which had been proven previously [27,28].To determine the regions of interest (ROIs) of task-evoked responses, we first computed the block-averaged motor response with local regressors (the standard) for each L-channel and subject.Then based on the block-averaged motor responses, we calculated the mean values from 5 to 20 seconds after task onset as the amplitudes of task-evoked responses.Finally, we selected two channels from each hemisphere that had maximum amplitudes as the ROIs for each subject and for each of the HbO and HbR signals.To evaluate the denoising performance of linear regression, within the identified ROIs, we calculated: (1) the Pearson's correlation coefficient (R) between the denoised real-time time series with local regressors (the standard) and the denoised time series with global regressors for each individual subject and (2) the mean square error (MSE) between them for each individual subject.The Fisher's transformation was then applied to convert R values to Z values: The mean values and the standard deviations of Z and MSE values were computed for each of the three types of global regressors: the Sym-Mid S-ch, Posterior S-ch, and Anterior S-ch regressors.Then the Z values were inversely transformed to R values.Finally, the differences between different types of global regressors were estimated by applying a twotailed t test on their Fisher-transformed R values and MSE values.

Spatial magnitude patterns of superficial interference
Figure 3 shows the power (N = 8) of the four physiological fluctuations in the anterior, middle, lateral and posterior regions, which was calculated based on the resting-state data in the first experiment.The power of very low-frequency fluctuations, Mayer waves and cardiac pulsations were of the same order, and all such fluctuations were much greater than the power of the respiratory waves.All four physiological fluctuations had relatively homogeneous magnitudes across the four scalp regions, and no spatial variation was seen at a significance level of p < 0.01.Furthermore, for all of the four physiological fluctuations, the power of HbO-related components was much greater than that of HbR-related components, indicating that these physiological fluctuations were primarily related to changes in arterial blood flow.Fig. 3. Power (mean ± standard error) of the four physiological fluctuations measured by the Schannels in multiple scalp regions: (A) power of HbO-related components, and (B) power of HbR-related components.In each graph, "A" represents the anterior region, "L" represents the lateral region, "M" represents the middle region, and "P" represents the posterior region.

Spatial similarity patterns of superficial interference
Figure 4 shows the spatial similarity patterns of superficial interference at both the local level and the global level, which was calculated based on the whole-frequency data in the first experiment.For HbO, the local-level similarities were very high in each of the anterior, middle, lateral and posterior regions (0.91 ± 0.05).The global-level similarities were high in the symmetrical pairs (0.76 ± 0.11), moderate in the ipsilateral (0.56 ± 0.15) and in the contralateral-asymmetrical pairs (0.56 ± 0.15).The t test at the group level (N = 8) showed that the local-level similarities were significantly higher than all three types of global-level similarities (p < 0.001, Bonferroni corrected).Furthermore, among the three types of globallevel similarities, the symmetrical pairs had significantly higher similarities than those of asymmetrical (both ipsilateral and contralateral-asymmetrical) pairs (p < 0.001, Bonferroni corrected).The difference between the ipsilateral and contralateral-asymmetrical similarities was insignificant (p > 0.05).For HbR, the local-level similarities were high at each of the anterior, middle, lateral and posterior regions (0.84 ± 0.09); the symmetrical similarities were moderate (0.54 ± 0.19), and the ipsilateral (0.34 ± 0.19) and contralateral-asymmetrical similarities (0.34 ± 0.19) were low.The local-level similarities were significantly higher than all three types of global-level similarities (p < 0.001, Bonferroni corrected).The difference between symmetrical similarities and asymmetrical (both ipsilateral and contralateralasymmetrical) similarities was also significant (p < 0.001, Bonferroni corrected).The difference between the ipsilateral and contralateral-asymmetrical similarities was insignificant (p > 0.05).Figure 5 shows the spatial similarity patterns of the four typical physiological fluctuations, which were calculated based on the band-pass filtered HbO data in the first experiment.All of these physiological fluctuations demonstrated decays of similarity from the local-level pairs to the global-level pairs, and from the globally symmetrical pairs to the globally asymmetrical (both ipsilateral and contralateral) pairs.All of these decays were statistically significant (p < 0.001, Bonferroni corrected).The difference between the ipsilateral pairs and the contralateral-asymmetrical pairs was statistically insignificant (p > 0.05).Specifically, the cardiac pulsations had very high similarity values at both the local level and the global level (local = 0.99 ± 0.01, symmetrical = 0.96 ± 0.03, ipsilateral = 0.94 ± 0.04, and contralateralasymmetrical = 0.94 ± 0.04).For HbR data, the spatial similarity patterns of the four physiological fluctuations are shown in Fig. 6.All of them demonstrated significant decays of similarity from the local-level pairs to the global-level pairs and from the globally symmetrical pairs to the globally asymmetrical pairs (p < 0.001, Bonferroni corrected), which were consistent with the results based on HbO data.Furthermore, the similarities calculated based on HbR data were lower than their counterparts based on HbO data.It was presumably because the HbR-related components in these physiological fluctuation had lower power (Fig. 3).As a result, the extracted HbR data were more vulnerable to measurement noise which subsequently reduced the similarity values.7.For the Mayer waves, respiration and cardiac pulsations, the similarities between the anterior region and all the other regions were lower than the similarities between any two non-anterior regions (the lateral, middle, and posterior regions).This difference was statistically significant for the respiration (p < 0.01, Bonferroni corrected) and cardiac pulsations (p < 0.001, Bonferroni corrected).For the Mayer waves, the similarity difference between the anterior region and the non-anterior regions also tended to be significant in HbO data.Overall, these results indicate that these three physiological fluctuations in the anterior region are relatively "isolated" from the other scalp regions.
The time lags of superficial interference among all of the local-level and global-level Schannel pairs were also investigated.For every S-channel pair, we recorded the absolute time lag (ranging from 0 to 50 sec) that showed the maximum cross-correlation coefficient.As summarized in Table 1, for the whole-frequency data, the absolute time lags were nearly zero among the local-level pairs, and gradually increased among the symmetrical pairs and the asymmetrical (both ipsilateral and contralateral-asymmetrical) pairs.Similarly, the absolute time lags for the low-frequency oscillations, Mayer waves and respiratory waves also demonstrated an overall increase from the local pairs to the symmetrical pairs and then to the asymmetrical pairs.The cardiac pulsations did not show remarkable time lags across different scalp regions, especially in HbO data.In Table 1, the absolute time lags estimated from HbR data consistently show greater mean values and standard deviations than their counterparts from HbO data.A close inspection showed that the HbR data were overall noisier than the HbO data due to their difference in power (referred to Fig. 3).Thus, the absolute time lags estimated from HbR data were more variable.

Denoising performance of linear regression
As a representative example, Fig. 8 shows the block-averaged motor responses in the second experiment for one subject, which were averaged across the L-channels within the identified ROI.A robust increase of HbO and decrease of HbR in the sensorimotor cortices are seen during the period of finger tapping.In particular, the block-averaged time courses with the local regressors (the standard) and the block-averaged time courses with the symmetricalmiddle regressors are very close to each other; both return to baseline after the task.In contrast, the block-averaged time courses with the bilateral posterior regressors and bilateral anterior regressors are inconsistent with the standard and do not return to baseline after the task.Figure 8 also shows the block-averaged motor responses by using the global mean of all S-channels as a regressor, which are close to those with the local regressors and symmetricalmiddle regressors and better than those with the bilateral posterior and anterior regressors.At the group level (N = 9), Fig. 9 shows the denoising performance of the four types of global regressors (symmetrical-middle, bilateral posterior, bilateral anterior and global mean) by comparing with the standard time series denoised by the local regressors.For HbO data, the symmetrical-middle regressors had the best denoising performance with highest R values and lowest MSE values; the bilateral anterior regressors had the worst denoising performance with lowest R values and highest MSE values.The denoising performance of the bilateral posterior regressors and global mean regressor was moderate.For HbR data, the symmetricalmiddle regressors also had the best denoising performance whereas the global mean regressor had the worst denoising performance.The denoising performance of the bilateral posterior regressors and bilateral anterior regressors was moderate.Furthermore, Fig. 9 also shows linear regression has overall worse denoising performance on HbO data.It might be because the physiological fluctuations are primarily related to changes in arterial blood blow (Fig. 3), which have greater influence on HbO data.

Discussion
In the current study, we investigated the similarity of superficial interference in multiregional fNIRS imaging in a whole-frequency domain (< 2 Hz) as well as in four specific frequency bands.In the whole-frequency domain, two major findings are summarized as follows: First, superficial interference showed very high similarities at the local level, whereas its similarities at the global level are significantly lower.This finding is consistent with the reports by Gagnon and his associates [27], which showed the correlation of superficial interference within a local area (< 1.5 cm) was greatest while the correlation of superficial interference between the motor area and the frontal area was much lower.Second, among the three types of S-channel pairs at the global level, the symmetrical pairs were found to have moderate to high similarities, which were significantly higher than those from the asymmetrical (both ipsilateral and contralateral-asymmetrical) pairs.This finding agrees with the previous reports about a symmetrical correlation pattern of fNIRS data [24] as well as blood pressure variance and heart rate variance in the same study.
In a further step, we investigated the similarity of superficial interference in four specific frequency bands that are related to the very low-frequency fluctuations, Mayer waves, respiration, and cardiac pulsations.All of these physiological fluctuations showed a gradual decay in similarity from the local S-channel pairs to the globally symmetrical pairs and then to the globally asymmetrical pairs.In particular, the cardiac pulsations showed very high similarities across all of the scalp regions in terms of HbO-related components.HbR-related components in the cardiac pulsations were not considered since the cardiac pulsations mainly induce changes in arterial blood flow, which is nearly 100% saturated.This finding indicates that the cardiac pulsations are globally homogeneous over the entire head, which also agrees with previous reports [15,27].
Interestingly, for the Mayer waves, respiratory waves and cardiac pulsations, we found that the similarities of the anterior region with all the other regions were significantly low.This phenomenon is not likely due to measurement artifacts because the power of these three physiological fluctuations in the anterior region is largely the same as in the other regions (Fig. 3).It implies that these three physiological fluctuations in the anterior region are relatively "isolated" from the other regions.The anatomical basis of this phenomenon needs to be further studied.A study by Kirilina et al. [29] has shown that task-evoked superficial artifacts in the prefrontal cortex are not homogeneously distributed, but are instead localized in the scalp draining veins.Their results are a helpful contribution to the study of the anatomical origin of superficial interference in the anterior region.Furthermore, although such a phenomenon is not found in the very low-frequency fluctuations, it does not mean the very low-frequency fluctuations are more homogeneous between the anterior region and the other regions.In fact, we noted from Fig. 7 that the cross-region similarities of the very lowfrequency fluctuations were lower than the other three physiological fluctuations.Therefore, the very low-frequency fluctuations are indeed more inhomogeneous in the human head.This study's observed spatial patterns of similarity pertaining to superficial interference provide an explanation as to why adjacent and remote S-channels lead to different denoising performance in regression of superficial interference, which has been noted but not explained in previous reports [27].In general, the hemodynamic activities in the brain are slow, and operate within a very low-frequency band (0.01-0.05 Hz) for a block-designed experiment and in the low-frequency band (0.05 -0.15 Hz) for an event-related experiment.Moreover the resting-state studies revealed a frequency range of functional connectivity from about 0.01 to 0.1 Hz [36][37][38].The results in this study demonstrate all three slow physiological fluctuations (very low-frequency fluctuations, Mayer waves and respiratory waves) with a gradual similarity decay from the local S-channel pairs to the globally symmetrical S-channel pairs and then to the globally asymmetrical (ipsilateral and contralateral-asymmetrical) Schannel pairs.Therefore, the different denoising performance of linear regression obtained by using different types of S-channels is reasonable and expected.This conclusion was further validated by the second experiment.
The results from this study also indicate that it is possible to simplify the configurations of S-channels in multiregional fNIRS imaging.As discussed in the beginning of this paper, the co-located configuration of S-channels and L-channels has proven to be optimal in regression of superficial interference.However, such a configuration is difficult to achieve in multiregional fNIRS imaging because it requires a large number of S-channels to be colocated with L-channels.Given the limited number of measurement channels that the current fNIRS systems have, reducing and simplifying the S-channels will always be an objective.A simplified fNIRS probe configuration will also reduce the preparation time and discomfort of subjects during experiments.According to the spatial similarity patterns of superficial interference identified in this study, superficial interference also has moderate to high similarity across the globally symmetrical S-channel pairs in addition to the local S-channels.Therefore, the globally symmetrical regions could share S-channels while maintaining an acceptable efficiency in regression of superficial interference.In this way, the number of Schannels can be reduced.
A few limitations in the current study should be noted.First, while a short source-detector separation (ideally ≤1 cm) for S-channels is always desired to ensure perfect specificity to the superficial layer [27], our S-channel separation of 1.5 cm was constrained by the optode size in the CW6 system.Nevertheless, the 1.5-cm separation of S-channels has been used in previous fNIRS studies to regress out superficial interference [16,17].Based on a finite element mesh (FEM) head model [39], we calculated the sensitivities of S-channels and Lchannels to the brain tissues.We found that the averaged sensitivity of S-channels to the brain was only 0.47% ± 0.14% when source-detector separation ranged from 1.4 to 1.6 cm.In contrast, the averaged sensitivity of L-channels to the brain was 3.55% ± 1.55% when sourcedetector separation ranged from 2.8 to 3.2 cm.Therefore, we are confident that the S-channels in the current study primarily measured superficial interference rather than brain activities.Furthermore, the very low sensitivities of S-channels to the brain indicate that linear regression of superficial interference using these S-channels will not reduce the amplitude of brain activities measured by the L-channels.A further demonstration on this issue is seen in the Appendix.Second, we extracted the typical physiological fluctuations from the scalprecorded fNIRS data according to their frequency bands.If we could record the peripheral physiological signals directly with auxiliary devices (such as pulse oximetry, respiratory belt and blood pressure monitor) and correlate them with the present fNIRS data from multiple scalp regions, it would enhance and expedite current comprehension of the spatiotemporal properties of these physiological fluctuations [24,40].This is a planned objective for future work.Lastly, the current study only measured superficial interference in selected scalp regions rather than the entire head, and only investigated the task-evoked activations at the sensorimotor cortex.In future research, the results from the current study will be replicated and validated by covering even larger scalp regions and by inducing other brain activations besides the sensorimotor cortex.over the sensorimotor cortices could definitely detect certain brain activations in the second experiment, it is plausible to think that linear regression with local S-channels might remove more true brain activations than the non-local S-channels and the global mean regressor.Because the results derived from the real experiment were lack of ground truth in terms of task-evoked brain activations, we tested this scenario in a simulative experiment.
In the simulative experiment, the real physiological noises were represented by the resting-state data in the second experiment, which was acquired prior to the motor tasks.The true brain activations were modeled by convolving a synthetic hemodynamic response function (HRF) with the stimulation function.In this way, a real controlled environment could be tested.Specifically, the synthetic HRF was generated as the sum of two Gamma functions.The amplitude ratio of the first Gamma function (peak time = 5 seconds) to the second Gamma function (peak time = 15 seconds) was 6:1, same as the canonical HRF in SPM8 (http://www.fil.ion.ucl.ac.uk/spm/).The amplitude of the synthetic HRF was set to be a 0.08-μM increase in HbO and a 0.023-μM decrease in HbR, which were estimated from the deconvoluted HRF in the motor tasks of the second experiment.The synthetic impulse responses were then convolved with an onset/offset time course of the motor tasks to yield the expected task-evoked hemodynamic responses in the brain (that is, the "ground truth").The onset/offset time course of tasks here was identical to the real motor task in the second experiment.Then, we selected four L-channels, two at the center of left sensorimotor cortex and the other two at the center of right sensorimotor cortex, as the ROIs.The simulated hemodynamic responses were added to the HbO and HbR data that were acquired from the four L-channels in the ROIs in the resting state.To simulate the small portion of brain responses detected by the local S-channels, we also added the simulated hemodynamic responses to four S-channels within the ROIs.Please note that, according to our calculation based on a FEM head model, the sensitivities of S-channels to the brain were about 1/7.5 of the sensitivities of L-channels to the brain.Therefore, the amplitudes of simulated hemodynamic responses added to the S-channels were 1/7.5 of those added to the L-channels.
For the simulated data series, linear regression was conducted on each L-channel in the ROIs by using every available S-channel as a regressor, same as it was carried out in the real experiment.The denoising performance of linear regression was evaluated with the R values and MSE values.At the group level (N = 9), Fig. 11 demonstrates that the local regressors did not reduce the amplitudes of recovered task-evoked brain response obviously.Figure 12 shows the results of R values and MSE values between the true and recovered task-evoked brain responses across all subjects (N = 9).The local regressors appear to have the best denoising performance, although their differences to the other regressors are statistically insignificant (p > 0.1).These results implies that superficial interfernce is the dominat component in the Schannels data series.Therefore, the linear regression operation mainly fitted and removed superficial interfernce from the L-channel data series.In summary, the results from the simulative experiment demonstrated that although the Schannels at 1.5-cm source-detector separation contained certain components from the brain, it did not lead to reduction of true brain activation from L-channels, and the denoising performance of linear regression with the local regressors was still the best.

Fig. 1 .
Fig. 1.Geometry of fNIRS probe and types of S-channel pairs in the first experiment: (A) Schematic arrangement of the sources, detectors and channels in the probe to record superficial interference at multiple regions over the head; (B) Illustration of the local-level S-channel pairs.(C) Illustration of the three types of global-level S-channel pairs: symmetrical pairs (sySch), ipsilateral pairs (ipS-ch), and contralateral-asymmetrical pairs (caS-ch).

Fig. 4 .
Fig. 4. Similarities of superficial interference in the whole-frequency domain (<2 Hz) measured by S-channels in multiple scalp regions.(A) Spatial distributions of S-channel pairs with similarities > 0.4 that were calculated based on HbO data.The colored lines represent the S-channel pairs in different ranges of similarity values.(B) Spatial distributions of S-channel pairs with similarities > 0.4 that were calculated based on HbR data.The colored lines represent the S-channel pairs in different ranges of similarity values.(C) Grand-averaged similarities of superficial interference at the local level and global level based on HbO data.(D) Grand-averaged similarities of superficial interference at the local level and global level based on HbR data.The error bars in (C) and (D) represent the standard deviations.

Fig. 5 .
Fig. 5. Similarities of superficial interference in four specific frequency bands based on HbO data: (A) the very low-frequency fluctuations between 0.01 and 0.05 Hz, (B) the Mayer waves between 0.07 and 0.15 Hz, (C) the respiratory waves between 0.2 and 0.4 Hz, and (D) the cardiac pulsations between 0.8 and 2 Hz.In each respective graph, the colored lines represent the S-channel pairs in different ranges of similarity values, and the error bars represent the standard deviations.

Fig. 6 .
Fig. 6.Similarities of superficial interference in four specific frequency bands based on HbR data: (A) the very low-frequency fluctuations between 0.01 and 0.05 Hz, (B) the Mayer waves between 0.07 and 0.15 Hz, (C) the respiratory waves between 0.2 and 0.4 Hz, and (D) the cardiac pulsations between 0.8 and 2 Hz.In each respective graph, the colored lines represent the S-channel pairs in different ranges of similarity values, and the error bars represent the standard deviations.

Figures 5
Figures 5 and 6 also indicate that the global-level similarities of the four physiological fluctuations across different scalp regions are quite different.Thus, we further evaluated the cross-region similarities of the four physiological fluctuations, which are shown in Fig.7.For the Mayer waves, respiration and cardiac pulsations, the similarities between the anterior region and all the other regions were lower than the similarities between any two non-anterior regions (the lateral, middle, and posterior regions).This difference was statistically significant for the respiration (p < 0.01, Bonferroni corrected) and cardiac pulsations (p < 0.001, Bonferroni corrected).For the Mayer waves, the similarity difference between the anterior region and the non-anterior regions also tended to be significant in HbO data.Overall, these results indicate that these three physiological fluctuations in the anterior region are relatively "isolated" from the other scalp regions.The time lags of superficial interference among all of the local-level and global-level Schannel pairs were also investigated.For every S-channel pair, we recorded the absolute time lag (ranging from 0 to 50 sec) that showed the maximum cross-correlation coefficient.As summarized in Table1, for the whole-frequency data, the absolute time lags were nearly zero among the local-level pairs, and gradually increased among the symmetrical pairs and the asymmetrical (both ipsilateral and contralateral-asymmetrical) pairs.Similarly, the absolute time lags for the low-frequency oscillations, Mayer waves and respiratory waves also demonstrated an overall increase from the local pairs to the symmetrical pairs and then to the

Fig. 7 .
Fig. 7. Cross-region similarities of superficial interference (the two paired S-channels were located in two different scalp regions).(A) and (C): The frequency-specific mean similarity matrices that were based on the HbO and HbR data, respectively.As shown in the legend, the upper triangular part of the similarity matrix is for the regions within the same hemisphere (ipsilateral pairs), and the lower triangular part is for the regions between two hemispheres (contralateral pairs).(B) and (D): Comparisons between the cross-anterior similarities (one of the two paired channels in similarity calculation was located in the anterior region) and notcross-anterior similarities (neither of the two paired channels in similarity calculation was located in the anterior region).The error bars in (C) and (D) represent the standard deviations.

Fig. 8 .
Fig. 8. Block-averaged motor responses from the L-channels within ROI for a representative subject.In each graph, the solid curves represent the block-averaged time courses, and the shaded areas indicate the standard errors across the seven blocks.Please note that the duration of finger tapping varied from 20 to 30 seconds among the seven blocks.

Fig. 9 .
Fig. 9. Denoising performance of linear regression in the second experiment: (A) R values between the denoised HbO time series using the four global regressors and the standard HbO time series using the local regressors.(B) MSE values between the denoised HbO time series using the four global regressors and the standard HbO time series.(C) and (D) R and MES results from the HbR data.The error bars in each graph represent the standard deviations.

Fig. 10 .
Fig. 10.Block-averaged motor responses from the four L-channels within the ROIs for a representative subject.In each graph, the red line represents the true task-evoked brain response, the non-red line represents the recovered task-evoked brain response after line regression, and the shaded areas indicate the standard errors of the recovered brain response across the seven blocks.

Figure 10
Figure10shows the block-averaged motor responses in the simulative experiment for a representative subject.They were averaged across the four L-channels within the ROIs.The recovered task-evoked brain responses with the local regressors and the symmetrical-middle regressors are very close to the ground truth.In contrast, the recovered task-evoked brain

Fig. 11 .
Fig. 11.Grand-averaged motor responses from the L-channels within ROIs for all the nine subjects.In each graph, the red line represents the true task-evoked brain response, the non-red line represents the recovered task-evoked brain response after line regression, and the shaded areas indicate the standard errors of the recovered brain response across the seven blocks.

Fig. 12 .
Fig. 12. Denoising performance of linear regression with different regressors in the simulative experiment: (A) and (C) R values between the true task-evoked brain responses and the recovered task-evoked brain responses; (C) and (D) MSE values between the true task-evoked brain responses and the recovered task-evoked brain responses.