Reduction of global interference of scalp-hemodynamics in functional near-infrared spectroscopy using short distance probes

Functional near-infrared spectroscopy (fNIRS) is used to measure cerebral activity because it is simple and portable. However, scalp-hemodynamics often contaminates fNIRS signals, leading to detection of cortical activity in regions that are actually inactive. Methods for removing these artifacts using standard source-detector distance channels (Long-channel) tend to over-estimate the artifacts, while methods using additional short source-detector distance channels (Short-channel) require numerous probes to cover broad cortical areas, which leads to a high cost and prolonged experimental time. Here, we propose a new method that effectively combines the existing techniques, preserving the accuracy of estimating cerebral activity and avoiding the disadvantages inherent when applying the techniques individually. Our new method accomplishes this by estimating a global scalp-hemodynamic component from a small number of Short-channels, and removing its influence from the Long-channels using a general linear model (GLM). To demonstrate the feasibility of this method, we collected fNIRS and functional magnetic resonance imaging (fMRI) measurements during a motor task. First, we measured changes in oxygenated hemoglobin concentration (∆Oxy-Hb) from 18 Short-channels placed over motor-related areas, and confirmed that the majority of scalp-hemodynamics was globally consistent and could be estimated from as few as four Short-channels using principal component analysis. We then measured ∆Oxy-Hb from 4 Short- and 43 Long-channels. The GLM identified cerebral activity comparable to that measured separately by fMRI, even when scalp-hemodynamics exhibited substantial task-related modulation. These results suggest that combining measurements from four Short-channels with a GLM provides robust estimation of cerebral activity at a low cost.


Introduction
Functional near-infrared spectroscopy (fNIRS) is a noninvasive functional neuroimaging technique that can measure concentration changes in oxygenated and deoxygenated hemoglobin (Δ Oxy-and Δ Deoxy-Hb) in the cerebral cortex. It has advantages of portability, fewer physical constraints on the participant, and simplicity of use. Therefore, although measurements are limited to the cortical surface, it has been adopted widely in clinical practices and daily life situations (Fujimoto et al., 2014;Hoshi, 2005;Kato et al., 2002;Obrig, 2014;Piper et al., 2014;Takeda et al., 2007;Vernieri et al., 2006).
However, undesirable artifacts such as head motion and scalphemodynamics often contaminate fNIRS signals, and obscure taskrelated cerebral-hemodynamics (for review, Scholkmann et al., 2014). In particular, scalp-hemodynamics, which is systemic changes of blood flows in the scalp layer, cannot be prevented experimentally because they are affected by systemic physiological changes resulting from activation of the autonomic nervous system or by changes in blood pressure accompanied by actions (Bauer et al., 2006;Lee et al., 2002;Scremin and Kenney, 2004). Indeed, both scalp-and cerebralhemodynamics increase in a task-related manner. This is especially true for ΔOxy-Hb, which is more widely used as an indicator of cerebral activity than ΔDeoxy-Hb because of its higher signal-to-noise ratio. For example, a majority of task-related changes in ΔOxy-Hb during a verbal fluency task was reported to originate from the scalp rather than the cortex (Takahashi et al., 2011). Furthermore, Minati et al. (2011) reported that rapid arm-raising movement during visual stimulus presentation generated transient increases in systemic blood pressure, and that Δ Oxy-Hb in the visual cortex was coupled with this change in blood pressure, rather than with visual stimulation. Given that scalp-and cerebral-hemodynamics in ΔOxy-Hb have similar temporal profiles, removing the scalp-hemodynamic artifacts by conventional temporal filtering or block averaging is difficult.
Assuming that changes in scalp-hemodynamics are more global than changes in cerebral-hemodynamics, several analytical techniques have been proposed that estimate scalp-hemodynamic artifacts from spatially uniform components of Δ Oxy-Hb that are measured by a standard source-detector distance of 30 mm (Long-channels). Using principal component analysis (PCA), Zhang et al. (2005) proposed an eigenvector-based spatial filter from data obtained during rest periods (baseline), which assumes that the effects of systemic hemodynamics is dominant in baseline data. This method has been further extended by applying Gaussian spatial filtering (Zhang et al., 2016). Furthermore, using independent component analysis (ICA), Kohno et al. (2007) extracted the most spatially uniform component of Δ Oxy-Hb and showed that it was highly correlated with scalp blood flow that was simultaneously measured by laser-Doppler tissue blood-flow. By removing these spatially uniform ΔOxy-Hb components, both methods identified task-related cerebral-hemodynamics in more spatially localized regions, suggesting that global scalp-hemodynamics is a major source of artifacts that decrease the signal-to-noise ratio in fNIRS measurements. Because the Δ Oxy-Hb recorded by Longchannels is a summation of scalp-and cerebral-hemodynamics, these techniques can lead to over-estimation of scalp-hemodynamic artifacts and underestimation of cerebral activity if the two spatially overlap or are highly correlated with each other. Therefore, independent measurement of scalp and cerebral hemodynamics-related Δ Oxy-Hb is preferable. Moreover, few studies have experimentally supported the assumed homogeneity of scalp-hemodynamics.
Recent studies have proposed removal of local scalp-hemodynamic artifacts using direct measurements from source-detector distances that are shorter than the standard Long-channels (these are the Shortchannels) (Funane et al., 2014;Gagnon et al., 2011Gagnon et al., , 2012aGagnon et al., , 2014Gregg et al., 2010;Saager et al., 2011;Yamada et al., 2009;Zhang et al., 2009Zhang et al., , 2011 (Fig. 1A). For example, Yamada et al. (2009) added a Short-channel detector with a 20-mm distance to each of four Longchannel probe pairs during a finger-tapping task, and subtracted the Short-channel signal from the corresponding Long-channel signal. They confirmed that the activation area that remained after artifact subtraction was comparable with that measured by functional magnetic resonance imaging (fMRI). Although this is a powerful and accurate technique, numerous probes are necessary to cover broad cortical areas because the same number of Long-and Short-channels is required. Dense and broad fNIRS probe arrangements are expensive, heavy, and time consuming, and therefore not practical or feasible, especially for clinical applications.
To take advantage of its simplicity of use and the ability to measure activity from broad cortical regions, a simple method that can remove fNIRS artifacts from broad measurement areas is required. To meet this requirement, here we first tested whether it is possible to estimate scalp-hemodynamic artifacts with a reduced number of Short-channels. Above-mentioned previous methods using Shortchannels used multiple Long-and Short-channel pairs because scalphemodynamics was considered to vary at different locations (channels) on the head. However, if scalp-hemodynamics is globally uniform as has been assumed in previous studies (Kohno et al., 2007;Zhang et al., 2005, Fig. 1. fNIRS probe arrangements. (A) Schematic illustration of depth sensitivity corresponding to source-detector distance. (B) Probe arrangement for the scalp-hemodynamics measurement in Experiment 1A. Eighteen Short-channels were arranged so that they covered the motor-related areas of both hemispheres. (C) Probe arrangement for evaluation in Experiment 2. Forty-three Long-channels and four Short-channels were arranged to cover the same areas as in (B). 2016), we can reduce the number of Short-channels necessary for estimating global scalp artifacts. Because the distribution of scalphemodynamics over broad measurement areas remains unclear, we first measured scalp-hemodynamics during a finger-tapping task using 18 Short-channels placed on bilateral motor-related areas, including primary sensorimotor, premotor, and supplementary motor areas (Experiment 1A) and during a verbal fluency task using eight Shortchannels on the prefrontal cortex (Experiment 1B). We assessed scalp-hemodynamic homogeneity and determined the number of Short-channels necessary for estimating artifacts.
Next, we tested a new method for estimating cerebral activity that combines minimal Short-channel measurements with a general linear model (GLM), a combination that has not yet been effectively employed toward broad fNIRS measurements. Results from Experiment 1A provided the minimal number of Short-channels needed to estimate global scalp-hemodynamic artifacts independently from cerebralhemodynamics in motor-related areas. The estimated global scalphemodynamic model was then incorporated into the GLM designmatrix together with the functional cerebral-hemodynamic models. Because simultaneous estimation of both scalp-and cerebralhemodynamics has been demonstrated to improve performance (Gagnon et al., 2011), our GLM with the scalp-hemodynamic model is expected to estimate cerebral activity more accurately by reducing the interference from scalp-hemodynamics. Although removal of scalphemodynamics is thought to be unnecessary or achievable through conventional means when scalp-hemodynamic artifacts and cerebral activity are independent, conventional methods cannot remove it when scalp-hemodynamics is highly correlated with cerebralhemodynamics.
Our study aimed to provide a practical approach that addresses scalp-hemodynamic artifacts by effectively combining already existing techniques, rather than developing a new and costly technique. We hypothesized that estimating global scalp artifacts using a few Short-channels would avoid the over-estimation occurs when using only Long-channels, and would be effective even when the scalphemodynamic artifact is correlated with cerebral-hemodynamics. We tested this hypothesis using data from fNIRS and fMRI experiments.

Proposed artifact reduction method
The proposed artifact reduction method consists of three steps: preprocessing, estimation of the global scalp-hemodynamic artifact, and removal of scalp-hemodynamics using GLM analysis. Specifically, in the second step, the common temporal pattern of scalphemodynamics was extracted from a small number of Short-channels using PCA. In the third step, the GLM with the scalp-hemodynamic (i.e., artifact) and cerebral-hemodynamic models in response to the task was applied to the Long-channels, yielding an estimation of cerebral activity free of global scalp-hemodynamic artifacts.

Preprocessing
According to standard practice, systemic physiological signals were filtered out from the fNIRS data if they occurred at frequencies that were higher or lower than once in a task cycle (45 s in this study) that did not correlate with the task cycle. The Δ Oxy-Hb signal from each channel was detrended using a discrete cosine transform algorithm with a cut-off period of 90 s (i.e. high-pass filtering with cut-off frequency of 1/90 Hz), which was twice the task cycle (see fNIRS data acquisition during the motor task section), and smoothed using a moving average with a temporal window of 3 s.

Estimation of global scalp-hemodynamics
Based on the results from Experiment 1A, the majority of scalphemodynamic artifacts were globally uniform (see below, Experiment 1A: distribution of scalp-hemodynamics in bilateral motor-related areas section). We therefore focused only on this global component and regarded it as the first principal component (1st PC) extracted from the Short-channel signals using PCA. We placed four Shortchannels on the bilateral frontal and parietal cortices at a probe distance of 15 mm (Fig. 1B) because this number was enough to precisely estimate the global component, according to the results from Experiment 1A (Fig. 3C). To prevent extracted PCs from being biased to a specific channel signal because of motion artifacts, signal-to-noise ratios, or different path lengths (Hoshi, 2005), the Δ Oxy-Hb signals were normalized for each channel before applying PCA so that the mean and standard deviation (SD) were 0 and 1, respectively.

Removal of scalp-hemodynamics using GLM
We applied a GLM to compute the relative contributions of scalpand cerebral-hemodynamics to each channel. We added the extracted hemodynamics model to the conventionally used GLM design matrix for fMRI (Friston et al., 2006) and fNIRS analyses (Koh et al., 2007;Plichta et al., 2006Plichta et al., , 2007Ye et al., 2009).
The GLM explains the response variable (fNIRS signals for each Longchannel) y in terms of a linear combination of the explanatory variables (design matrix) X: Here, β is an unknown weight parameter vector and ε is an error term that is independent and normally distributed with a mean of zero and variance σ 2 , i.e., ε~N(0, σ 2 ).
In the proposed model, design matrix X is constructed with five components: the cerebral-hemodynamic model to be predicted, its temporal and dispersion derivatives, a constant, and the global scalphemodynamic model (see Estimation of global scalp-hemodynamics section). Here, the cerebral-hemodynamic model is calculated by convoluting the task function and the hemodynamic response function (Friston et al., 2006). The temporal and dispersion derivatives can model small differences in the latency and duration of the peak responses, respectively (Friston et al., 1998).
The relative contributions of cerebral-and global scalphemodynamics to each Long-channel signal (i.e. estimated weight parameter β) are calculated using the least square method. We can evaluate the cerebral activity for each channel in each sample (see below) by calculating sample-wise t-statistics for β corresponding to the cerebral-hemodynamic model (Friston et al., 2006). Note that the t-statistic represents how much larger the estimated cerebralhemodynamics value is than the error, and is calculated at a single sample level in order to independently evaluate the estimation for each sample. The temporal correlation that exists in the residual error term was corrected using the precoloring method (Worsley and Friston, 1995). Preprocessing of the global scalp-hemodynamics model for the precoloring was performed before its extraction from the Short-channels. We can remove the global scalp-hemodynamics from each Long-channel signal by subtracting the global scalphemodynamic model multiplied by the β for the global scalphemodynamic model for the corresponding channel.
fNIRS data acquisition during the motor task In the motor task (Experiments 1A and 2), fNIRS measurements were performed using a multichannel continuous-wave optical imaging system (16 sources and 16 detectors) with wavelengths at 780, 805, and 830 nm (FOIRE-3000; Shimadzu Corp., Kyoto, Japan). Because we assume that the proposed method will be useful in clinical practices such as rehabilitation, we used 32 probes to cover the motor-related areas of both hemispheres (bilateral primary sensorimotor cortices, dorsal premotor cortices, and supplementary motor areas; Fig. 1B and C for Experiments 1A and 2, respectively). Measurements were performed during repetitive sets of motor tasks following a simple block design: six rest-task-rest (15 s each) blocks, with a 100 ms (Experiment 1A) or 130 ms (Experiment 2) sampling period. The start of the task period was indicated by a single clicking sound and the end by double clicks. Participants sat in a comfortable reclining armchair with both hands resting naturally on their knees, and were presented with a fixation point approximately 1 m in front of their faces. They were asked not to move their bodies during the rest period, and to repetitively tap their right index finger (brisk extension and flexion) or grasp a ball with either their left or right hand as fast as possible during the task period. After fNIRS measurement, each fNIRS probe position was recorded with a stylus marker (FASTRAK; Polhemus, Colchester, VT, USA).
We assumed that data obtained from the left-and right-hand tasks were independent, even when measured from the same participant, as in the method adopted in our previous study (Nambu et al., 2009). Thus, we analyzed each as a separate sample.
Although both ΔOxyand ΔDeoxy-Hb were calculated according to the modified Beer-Lambert Law (Delpy et al., 1988), we only analyzed ΔOxy-Hb because detecting ΔDeoxy-Hb with the wavelength pairs of our fNIRS system is difficult, and produces low signal-to-noise ratios (Hoshi, 2005;Mihara et al., 2012;Sato et al., 2004;Uludag et al., 2004) (see Discussion and Supplementary material).
Experiment 1A: scalp-hemodynamics measurement by multiple Shortchannels in the motor task First, we verified that the major component of the scalphemodynamics was consistent over a broad range of cortical areas, as has been assumed in previous studies (Kohno et al., 2007;Zhang et al., 2005Zhang et al., , 2016, and that it could be extracted from a few Short-channels. Thirteen right-handed healthy volunteers (aged 22-34 years, 12 men and one woman) participated in this experiment. Participants gave written informed consent for the experimental procedures, which were approved by the ATR Human Subject Review Committee and institutional review board of Nagaoka University of Technology.
To measure the scalp-hemodynamics in motor-related areas of both hemispheres during a hand movement task, 18 Short-channels were arranged so that the center of the probes corresponded to the Cz position of the International 10-20 system (Fig. 1B). The experimental protocol is given above in the fNIRS data acquisition during the motor task section. Three participants performed a tapping task with their right index finger (three samples), and the other 10 participants performed a grasping task twice, once with each hand (20 samples). Two samples were excluded from analysis because motion artifacts were found to have contaminated the raw fNIRS signals. Thus, 21 samples (8 for left-hand movements and 13 for right-hand movements) were analyzed.
After preprocessing, correlation coefficients of Δ Oxy-Hb on each Short-channel pair were calculated in each sample. Using PCA, the 1st PC was extracted from the 18 channels as the major component of the scalp-hemodynamics, and its contribution ratio was calculated for each sample. The contribution ratio is a measure of how much the PCs explain the variance, which is also referred to as the variance accounted for (VAF).
To examine the number of Short-channels necessary to reliably extract the global scalp-hemodynamics, we computed a correlation index (CI) as follows: Corr PC1 X short Here, A n denotes arbitrary channel combinations of n channels out of N ch Short-channels (e.g., A 4 = {1, 2, 3, 4}), X i Short indicates the Shortchannel signal obtained from channel i, PC1(X {A n } Short ) is the 1st PC extracted from the channel combination A n , and Corr(PC1,X i Short ) is the correlation coefficient between the extracted 1st PC and Shortchannel signal X i Short . In each sample, the CI was computed for all the possible combinations out of the 18 Short-channels and its median value for each n channel combination was obtained. Then, the average of median CI for each n was calculated across all the samples.

Experiment 1B: scalp-hemodynamics measurement in the verbal fluency task
To investigate the homogeneity of scalp-hemodynamics during a non-motor task, we conducted an experiment with a verbal fluency task and measured scalp-hemodynamics in the prefrontal cortex using Short-channels. We adopted the protocol used in Takahashi et al. (2011). The fNIRS probes for measurements of eight Short-channels were placed on the forehead (for details, see Supplementary material 2). Ten participants from Experiment 1A participated in this additional experiment, but four were excluded from the analysis because the data were contaminated by motion artifacts (checked by visual inspection). Preprocessing was the same as in Experiment 1A, except that detrending (high-pass filtering) was set to a cut-off period of 320 s (instead of 90 s) because the task cycle was 160 s (instead of 45 s).

Experiment 2: experimental evaluation of the performance of the proposed method Samples
We investigated the performance of the ShortPCA GLM using measurements from Long-and Short-channels. Sixteen right-handed volunteers, 15 healthy participants (aged 22-67 years, seven men, eight women; seven had also participated in Experiment 1A) and one right-handed male stroke patient (aged 60, 6 years after stroke), participated in this experiment. The patient had an infarction in the right corona radiata and showed mild left hemiparesis (Brunnstrom stage V for hands/fingers: voluntary mass finger extension with a variable range of motion; Brunnstrom, 1966). All participants gave written informed consent, and the experiment was approved by the ATR Human Subject Review Committee, institutional review board of Nagaoka University of Technology, and the ethical committee of Tokyo Bay Rehabilitation Hospital. fNIRS signals were recorded by the probe arrangement described above (Fig. 1C). Two healthy participants performed a tapping task with their right index fingers, and the others performed a grasping task twice, once with each hand, as described in the fNIRS data acquisition during the motor task section. Patient data obtained for the unaffected right hand were excluded from analysis because of motion artifacts that contaminated the raw fNIRS signals. Two healthy participants who performed the grasping task were also excluded because they were not relaxed during the fMRI task and unusual bilateral fMRI activation was observed (see Supplementary material 1 for details). Thus, 25 samples (12 for left-hand movements and 13 for right-hand movements) were analyzed.

Evaluation of proposed method
First, we calculated a variance inflation factor (VIF) (Kutner et al., 2004) to test multicollinearity in the GLM because VIF higher than 10 usually causes significant problems for estimating the parameters in GLM. Next, the performance of our new method was compared with those of three conventional methods: RAW, MS-ICA (ICA algorithm proposed by Molegedey and Schuster; Kohno et al., 2007), and RestEV (an eigenvector-based spatial filtering method using the rest; Zhang et al., 2005). The RAW method directly applied the Standard GLM to the preprocessed fNIRS signals. The MS-ICA and RestEV methods first estimated global artifacts using only Long-channels that covered wide cortical areas, and then removed the estimated artifacts from the fNIRS signals (see MS-ICA and RestEV sections for definitions and more detail). Cerebral-hemodynamics was then estimated by applying the Standard GLM.
To evaluate results of the GLM analysis, we first separated the fNIRS samples based on the degree of correlation between cerebraland scalp-hemodynamics. Samples in which the global scalphemodynamic model did not significantly correlate with the cerebralhemodynamic model (i.e., the influence of scalp-hemodynamics on the estimation of cerebral-hemodynamics was low) were assigned to the cerebral-scalp uncorrelated group. Those in which the global scalp-hemodynamic model significantly correlated with the cerebralhemodynamic model were assigned to the cerebral-scalp correlated group. The significance of the correlation was assessed by a permutation test. We computed the null distribution of the correlation as follows: for each sample, we first generated 100 sets of a "task-onset randomized" cerebral-hemodynamic model, which was calculated by convolving the hemodynamic response function and the boxcar task function that randomly changed task-onset. We then calculated the correlation coefficients between the global scalp-hemodynamic model of each sample and the randomized cerebral-hemodynamic model. Finally, we obtained the correlation threshold as the 95th percentile (one-tailed) of the null distribution for the correlation coefficients of 2500 samples (25 samples × 100).
To evaluate the goodness-of-fit of the GLM for each group, an adjusted coefficient of determination (adjusted R 2 ) was calculated and averaged over all Long-channels for each of RAW, MS-ICA, RestEV, and ShortPCA methods. For comparison, a two-tailed paired t-test was applied to assess differences in the averaged adjusted R 2 between ShortPCA and each of the other three methods (with Bonferroni correction, n = 3, significant level α = 0.017), respectively.
To confirm whether cerebral-hemodynamics estimated from the fNIRS signals accurately reflected cerebral activity, we performed an fMRI experiment with a similar task and compared fNIRS-estimated cerebral activity with that estimated by fMRI (see fMRI experiment section).

MS-ICA
Kohno et al. (2007) removed a component that had the highest coefficient of spatial uniformity among the independent components separated by an algorithm proposed by Molegedey and Schuster (referred to as MS-ICA; Molgedey and Schuster, 1994), and considered it the global scalp-hemodynamics component. Their artifact removal algorithm was implemented into the analytical software in our fNIRS system, and we applied it to the Δ Oxy-Hb measured from the 43 Long-channels. After removing the artifact, the signal was processed using the same preprocessing procedures (detrending and smoothing) as our new method, the Standard GLM was applied, and the results were compared with those obtained from the ShortPCA GLM.
RestEV Zhang et al. (2005) proposed an eigenvector-based spatial filtering method using the rest (baseline) period (referred to as RestEV). This method removes the first r spatial eigenvectors calculated from baseline data by PCA (for details, see Zhang et al., 2005). Here, the spatial filtering was applied to the preprocessed ΔOxy-Hb from the Long-channels, and then the Standard GLM was applied. The first rest period (15 s) was used to determine the eigenvector-based spatial filter. Zhang et al. (2005) determined the number of components r to be removed based on the spatial distribution of the eigenvector, although they did not give any clear criterion. We report the results of r = 1, as we confirmed that spatial filtering with a different r value (from 1 to 3) showed similar results (data not shown).

fMRI experiment
To confirm whether cerebral-hemodynamics estimated from fNIRS signals accurately reflected cerebral activity, we performed an fMRI experiment with a similar task and compared cerebral activity estimated by ΔOxy-Hb signals with that estimated by fMRI blood-oxygenation level dependent (BOLD) signals. This is reasonable given that ΔOxy-Hb signals are temporally and spatially similar to BOLD (Cui et al., 2011;Fig. 2. Results of the scalp-hemodynamics measurement (Experiment 1A). (A) ΔOxy-Hb from all 18 Short-channels (broken lines) and the 1st PC (thick line) for a representative sample. The signal of each channel was normalized. Gray shaded areas indicate the task period. (B) Waveform of block-averaged 1st PC, ensemble averaged across samples (thick line; thin lines indicate SD) and the cerebral-hemodynamic model (broken line). Each waveform was normalized so that the peak-to-peak amplitude was 1. (C) Average median values for the correlation indexes across samples for each number of Short-channels used to extract 1st PC. The error bar indicates the SD over samples. Okamoto et al., 2004b;Strangman et al., 2002;Toronov et al., 2007). We evaluated the estimation accuracy for each method using estimation error and signal detection measures, assuming that fMRI correctly detects cerebral activity. Thus, effectiveness of the proposed method was examined by comparing the estimation accuracy of the ShortPCA with the three conventional methods.
All participants of the fNIRS experiment also participated in the fMRI experiment, and performed the same motor task. T1-weighted structural images and functional T2*-weighted echo-planar images were recorded using a 1.5 or 3.0 T MRI scanner. Details for the fMRI experiment and parameters are described in Supplementary material 1.
Functional images were analyzed using SPM8 software (Wellcome Trust Centre for Neuroimaging, University College London, London, UK; http://www.fil.ion.ucl.ac.uk/spm/) to obtain parametric cerebral activation maps related to the task for each sample. First, images were preprocessed (realignment to correct head motion, co-registration of the functional images to each participants' structural image, and smoothing with a 9 mm full-width at half maximum Gaussian filter). Then, the Standard GLM was applied to the smoothed functional images to calculate the voxel t-values for the difference between the task and rest periods.
To compare fNIRS and fMRI activation maps for each sample, the fMRI t-value for the cerebral-hemodynamic model corresponding to each fNIRS channel location was calculated based on a previous report (Singh et al., 2005; details are described in Supplementary material 1). To normalize the difference in the degrees of freedom between fNIRS and fMRI GLMs, each t-value was divided by a t-value of the corresponding degree of freedom at a significant level with Bonferroni correction that was divided by the number of Long-channels (p = 0.05, n = 43, significant level α = 0.0012). On the assumption that the fMRI correctly reflects cerebral activity, we compared the t-values given by fNIRS with those given by fMRI.

Evaluation of estimation accuracy by comparing fMRI and fNIRS
The estimation accuracy for fNIRS activity for cerebral-scalp uncorrelated and correlated groups were evaluated by an estimation error and the signal detection measures computed using the corrected and aligned t-values described above. Note that these evaluations were conducted using sample-wise t-values (not group-averaged) because our purpose was to compare fMRI and fNIRS results for each sample. The estimation error was the root-mean-square difference between the fNIRS and fMRI t-values over all Long-channels. To compare the estimation error of the ShortPCA GLM with those of the other three methods (RAW, MS-ICA, and RestEV), we performed a two-tailed paired t-test (with Bonferroni correction, n = 3, significant level α = 0.017). As fNIRS analysis tends to estimate false positive cerebral activity, we examined the signal detection property of the ShortPCA GLM compared with the other methods. We calculated sensitivity, specificity, and the geometric mean of the sensitivity and specificity (G-mean) as the signal detection measures: where TP, TN, FP, and FN are the number of true positives, true negatives, false positives, and false negatives, respectively. Here, the true positive counts the number of channels in which both fNIRS and fMRI showed a significant increase in cerebral-hemodynamics, while the true negative counts the number of channels in which neither signal showed a significant increase (p = 0.05 with Bonferroni correction, n = 43, significant level α = 0.0012). The false positive counts the number of channels in which a significant increase in cerebralhemodynamics was detected by fNIRS but not by fMRI, while the false negative counts the number of channels in which the reverse was true. Sensitivity and specificity indicate the degrees to which actual activation and lack of activation are correctly detected, respectively. Gmean combines sensitivity and specificity, and is an appropriate measure when the number of actually activated and inactivated channels are imbalanced. As the number of actually activated channels was Estimated scalp-and cerebral-hemodynamics are represented by dotted and broken curves, respectively. The signal of each channel was normalized (z-score) by the mean and SD of the rest period before each task block.
limited, measures for each method were calculated using the pooled data in each group.

Results
Experiment 1A: distribution of scalp-hemodynamics in bilateral motorrelated areas Fig. 2A shows the ΔOxy-Hb for all 18 Short-channels (broken lines) and the 1st PC (thick line) for a representative sample of right-hand movement. The temporal changes in Δ Oxy-Hb were similar in each of the 18 Short-channels. The other samples showed a similar tendency. We calculated average correlation coefficients between all Short-channel pairs for each sample, and found that the mean of all 21 samples was 0.838 (SD 0.073). The major component of the scalphemodynamics was extracted by applying PCA to these data. The VAF of the 1st PC among all the samples was 0.850 (SD 0.067), and the mean value of the average correlation coefficients between the 1st PC and all Short-channel signals was 0.920 (SD 0.038). These results indicate that the 1st PC can explain the majority of the component contained in ΔOxy-Hb obtained from the 18 Short-channel signals.
To investigate the characteristics of task-related changes in global scalp-hemodynamics, the block-averaged 1st PC was ensemble averaged over all the samples (Fig. 2B, thick line). Compared with the cerebral-hemodynamic model (broken line), the 1st PC increased quickly after the onset of the task, reached a peak approximately 5 s earlier than the cerebral-hemodynamic model, and then decreased gradually after the peak. Fig. 2C shows the correlation index (CI, Eq. (2)). The average of median CI across samples increased exponentially as the number of Shortchannels used in the PCA increased. Importantly, the correlation index calculated for all 18 channels was similar to those for only 2, 3, and 4 channels. The differences were 0.031 (SD 0.016), 0.020 (SD 0.010), and 0.014 (SD 0.007), respectively. We also confirmed that a similar trend was produced using the ensemble average of ΔOxy-Hb over the selected Short-channels, instead of using their 1st PC.

Experiment 1B: distribution of scalp-hemodynamics in the verbal fluency task
Consistent with the results from Experiment 1A (Fig. 2), very similar temporal patterns in Δ Oxy-Hb were observed for all measured Short-channels (Suppl. Fig. 2A in Supplementary material 2). The average correlation coefficient across all channel pairs was 0.838 (SD 0.095) and the contribution ratio of the first principle component was 0.853 (SD 0.073). The mean value of the average correlation coefficients between the 1st PC and all Short-channels was 0.916 (SD 0.052). These values were similar to those obtained in Experiment 1A. The correlation index also showed a tendency similar to that observed in Experiment 1A (Suppl. Fig. 2C). The differences from the correlation index calculated for all eight channels were also relatively small for indexes calculated for two (0.026, SD 0.013), three (0.014, SD 0.006), and four channels (0.009, SD 0.004). These results suggest that scalp-hemodynamics in the prefrontal cortex during the verbal fluency task were as homogenous as those in the motor-related areas during the hand-movement tasks.

Evaluation of the model in real fNIRS signals
We evaluated performance of the ShortPCA GLM using experimental data, and compared the results with those obtained from three other conventional methods (RAW, MS-ICA, and RestEV). Fig. 3 shows the temporal waveform for the global scalp-hemodynamic model (Fig. 3A) and estimation results (Fig. 3B) using ShortPCA for a representative sample (left-hand, cerebral-scalp correlated group). The global scalp-hemodynamic model extracted from the four Short-channels showed a task-related temporal characteristic (Fig. 3A) such that hemodynamics increased during the task period and decreased during the rest period. The mean contribution ratio of the global scalphemodynamics component (1st PC) was 0.853 (SD 0.082). Among all samples, the maximum VIF we observed was 5.96, indicating that multicollinearity was not an issue for these samples. Fig. 3B shows an example of the estimated cerebral-and global scalp-hemodynamics from Ch 21 and Ch 27, placed on the left (ipsilateral) and right (contralateral) primary sensorimotor cortices, respectively. The majority of Δ Oxy-Hb (solid line) in Ch 21 was explained solely by global scalp-hemodynamics (Fig. 3B left). Conversely, Ch 27 contained both cerebral-and global scalp-hemodynamic components.
We next compared the results between cerebral-scalp uncorrelated and correlated groups. Using a correlation threshold of 0.314 (determined by permutation test), seven samples were assigned to the uncorrelated group while the remaining 18 samples were assigned to the correlated group.
To compare the fitting of each GLM, an adjusted R 2 was calculated for each of the Long-channels and averaged over the 43 Long-channels for each sample (Fig. 4). The adjusted R 2 for the proposed ShortPCA GLM was significantly higher than that for the other methods in the uncorre-  Fig. 4B). These results suggest that ShortPCA GLM is the most appropriate for fitting the Δ Oxy-Hb during movements in a block design. Fig. 5 shows the cerebral activity t-maps estimated by fNIRS (ΔOxy-Hb) and fMRI (spatially down-sampled to the locations of the fNIRS Long-channels) for group-level analysis of a right-hand sample assigned to the cerebral-scalp correlated group (Fig. 5A, C), and for a patient's affected left-hand sample, also assigned to the cerebral-scalp correlated group (Fig. 5B, D). The t-values of the estimated cerebralhemodynamics at the 43 Long-channels were mapped on the brain surface created by the structural MR image using Fusion software (Shimadzu Corp.;Okamoto et al., 2004a;Takeuchi et al., 2009) based on the fNIRS probe position recorded by a stylus marker.

Estimation accuracy in comparison to fMRI
For RAW, we found significant activity in all Long-channels. However, fMRI results demonstrated that brain activity was localized to the contralateral sensorimotor area, suggesting that false positive activity was estimated because of task-related changes in scalp-hemodynamics. Thus, without accurate removal of scalphemodynamics, cerebral activity can be over-estimated. Results applying scalp-hemodynamics removal showed improvements in the estimation of cerebral activity. Importantly, ShortPCA seemed to produce an activity map more similar to the fMRI result at the single sample level (Fig. 5B, D). In the representative samples assigned to the cerebral-scalp correlated group, when scalphemodynamics were estimated using MS-ICA or RestEV, and removed during preprocessing, the estimated cerebral activity was distributed in some unrelated areas. Additionally, a large negative value was estimated in some channels. In contrast, ShortPCA accurately estimated the expected cerebral activity, localizing it in the contralateral primary sensorimotor cortex. When we compared the fMRI t-maps with those from each fNIRS method, the ShortPCA spatial map showed similar activation patterns to those obtained by fMRI.
To quantify these results for each sample level, we evaluated the estimation accuracy for the cerebral-scalp uncorrelated and correlated groups. For the uncorrelated group, ShortPCA, estimation error did not significantly differ from the other methods (significant level α =   Fig. 6A). The signal detection measures were also similar (Fig. 6C). Conversely, estimation error for the correlated group (Fig. 6B), was significantly lower for ShortPCA than for the other methods (significant level α = 0.017 after Bonferroni correction; for vs RAW, paired t[17] = − 2.84, p = 0.0112; for vs MS-ICA, paired t[17] = −2.78, p = 0.0129; for vs RestEV, paired t[17] = −3.16, p = 0.00567). When we looked at the signal detection measures for the correlated group, RAW had the highest sensitivity, but the lowest specificity (Fig. 6D). The other three methods drastically improved the specificity compared with RAW. In particular, ShortPCA had the highest sensitivity and highest specificity among the three methods, and the highest G-mean among all four.

Discussion
In this study, we propose a new method (ShortPCA) for removing task-related scalp-hemodynamic artifacts that cannot be filtered out by conventional fNIRS analysis. Our method extracts the global scalphemodynamic artifact from four Short-channels, and then uses a GLM to simultaneously remove this artifact and estimate cerebral activity. Our method proved to be successful using fNIRS and fMRI experimental data. Among several alternative models, ShortPCA accurately fitted ΔOxy-Hb data obtained from fNIRS, and produced an estimated cerebral activity pattern that was the most similar to that observed by fMRI.
Our method is an effective combination of two previously developed techniques that have been independently studied. Studies for extracting global components have used only Long-channels, and no study has tried to incorporate Short-channel data toward this goal. Similarly, studies that directly measured scalp-hemodynamics using Shortchannels focused only on the accurate estimation of local artifacts, and did not attempt to estimate global components. By taking advantage of these two techniques, we were able to overcome their individual drawbacks (over-estimation of global scalp artifacts and measurement complexity).

Homogeneity of scalp-hemodynamics
Scalp-hemodynamics have been considered to be inhomogeneous in several studies using multiple Short-channels (Gagnon et al., 2011(Gagnon et al., , 2012aGregg et al., 2010;Saager et al., 2011;Yamada et al., 2009;Zhang et al., 2009Zhang et al., , 2011. Recently, Gagnon et al. (2012a) showed that the initial baseline correlation of fNIRS time courses between Shortchannels decreased with the increase in relative distance between the two channels, suggesting localization of scalp-hemodynamics during rest period. Furthermore, Kirilina et al. (2012) used fNIRS and fMRI and reported that scalp-hemodynamics were localized in the scalp veins, indicating a locally regulated physiological process in the scalp.
Other reports advocate homogeneity of scalp-hemodynamics (Funane et al., 2014;Kohno et al., 2007;Zhang et al., 2005Zhang et al., , 2016, which is physiologically supported by the reflection of systemic physiological changes in heart rate, respiration, and arterial blood pressures in fNIRS signals (Kirilina et al., 2013;Obrig et al., 2000). Our results regarding scalp-hemodynamics (Experiment 1A, Fig. 2A) revealed that ΔOxy-Hb measured by the Short-channels was uniform over the measured areas, and the majority of this temporal pattern was explained by the 1st PC. Similar spatially homogeneous characteristics of scalphemodynamics were observed by an additional experiment with Short-channel measurements of prefrontal cortex during a verbal fluency task (Experiment 1B). To our knowledge, this is the first experimental result that quantitatively supports homogeneity of scalp-hemodynamics over broad measurement areas. Although this result is inconsistent with the literature showing it to be inhomogeneous, this discrepancy is likely because of differences in experimental settings such as the measurement hardware (and the source-detector separation) rather than differences in the methods used to evaluate homogeneity. Indeed, we computed the initial baseline correlation between two Short-channels in the same way as Gagnon et al. (2012a) for all the samples, and verified that the correlation remained high even though the relative distance between the two channels increased. To clarify any limitations of our method, further studies should investigate the range of experimental conditions in which scalphemodynamics show homogeneity.

Required number of Short-channels
Based on the experimental evidence that the majority of scalphemodynamic artifacts are globally uniform, we can use fewer Short-channels to extract the scalp-hemodynamic artifact than past approaches have done. Despite the uniformity of the scalphemodynamics, we consider that more than one channel is required for the following reasons. First, each Short-channel signal may be contaminated with local noise which has similar frequency characteristics to the task cycle and therefore cannot be filtered out (e.g., local motion of the probe). Second, some channels occasionally fail to work properly because of poor contact between the probe and the scalp. Thus, separate location of Short-channels is generally better for applying our method to avoid local noises including motion artifacts.
The more Short-channels we used for extracting the global scalphemodynamics, the higher the correlation indexes were (Fig. 2C). Nevertheless, they were very high even with a few Short-channels. When the number of available probes is limited, arranging many Short-channels to improve the artifact estimation reduces the number of Long-channel probes and prevents coverage of broad cortical areas. Therefore, the minimum number of Short-channels is preferable. In Experiment 2, we assigned four Short-channels with the arrangement shown in Fig. 1C so that the motor-related areas are covered with Long-channels using 16 probe pairs. To examine the efficacy of this specific arrangement in estimating global scalp-hemodynamics, we compared the correlation index of this specific four Short-channel combination with that of all 18 Short-channels using the data from Experiment 1A. The difference in correlation indexes averaged across 21 samples was 0.017 (SD 0.013), suggesting that this arrangement was adequate for estimating cerebral-hemodynamics around motorrelated areas in our experimental setting. Location and combination of Short-channels should not be crucial because the difference between the highest and the lowest CI among all possible four Short-channel combination pairs was very small (around 0.03). However, the optimal number and location of Short-channels may vary depending on experimental conditions and the homogeneity of scalp-hemodynamics (Haeussinger et al., 2014;Kirilina et al., 2012), and care should be taken to evaluate them accordingly.

Comparison of performance evaluation using actual fNIRS signals
We compared the proposed ShortPCA method with three conventional methods (RAW, MS-ICA, and RestEV) by looking at the adjusted R 2 values (Fig. 4), t-value estimation errors (Fig. 6A, C), and the ability to detect cerebral activation (i.e. signal detection measures) (Fig. 6B, D). We found that ShortPCA performed best among all the methods when the scalp-hemodynamics significantly correlated with the cerebral-hemodynamics.
Scalp-hemodynamics often increases in a task-related manner, and over-estimation of the artifact and under-estimation of cerebral activity is an issue if only Long-channels are used to reduce the task-related artifact. Consistent with this idea, estimation accuracy degraded for conventional methods when applied to the cerebral-scalp correlated group. In the RAW method, large false-positive cerebral activity was observed in many channels as shown in Fig. 5, resulting in very low specificity (Fig. 6D). Although MS-ICA and RestEV were able to remove the scalphemodynamic artifact and improve specificity, they had lower adjusted R 2 in the GLM (Fig. 4) for both cerebral-scalp uncorrelated and correlated groups. These poorer fittings likely resulted from only using Longchannels, which led to misestimating the global scalp-hemodynamic components from fNIRS signals. As demonstrated by a simulation (Supplementary material 3), global scalp-hemodynamic models extracted from only the Long-channels were contaminated by cerebral-hemodynamics and caused the global scalp-hemodynamic model contribution to the GLM to be over-estimated. To support this, we found that G-mean was lower, owing to the lower sensitivity (representing false negatives), when we performed GLM analysis using only Long-channels (Supplementary material 8). This tendency was evident for the cerebral-scalp correlated groups. Thus, our results indicate that by using four Short-channels, robust estimation of cerebral activity can be achieved regardless of whether the influence of scalp-hemodynamics is significant.

Methodological considerations and limitations of the proposed method
To evaluate the accuracy of the estimated cerebral-hemodynamics, fMRI BOLD signals obtained during the same motor task were used as a standard for correct cerebral-hemodynamics. Because we did not measure fNIRS and fMRI simultaneously (owing to limitations of our fNIRS system), the cerebral activity measured in the two experiments was not identical. We cannot deny the possibility that the small differences that we observed between the two measurements were derived from differences in behavior, physical responses, or mental states. Additionally, although we compared t-values derived from Δ Oxy-Hb to those from fMRI-BOLD signals (as has been done previously: Cui et al., 2011;Okamoto et al., 2004b;Strangman et al., 2002;Toronov et al., 2007), which fNIRS signal correlate best with the BOLD signal remains an open question (Gagnon et al., 2012b;Huppert et al., 2006;Kleinschmidt et al., 1996;Mehagnoul-Schipper et al., 2002;Seiyama et al., 2004;Steinbrink et al., 2006;Toronov et al., 2001). Several studies indicate that during motor tasks, ΔOxy-Hb in frontal or ipsilateral areas is higher when measured by fNIRS than by fMRI (Cui et al., 2011;Okamoto et al., 2004b;Takeda et al., 2014). Thus, differences in detectability of cerebral activity could account for this aspect of our results. However, considering the robust and reproducible hemodynamic responses evoked by simple motor tasks (Yousry et al., 1997), and the spatial resolution of fNIRS that is not high enough to detect slight differences in protocols, we believe that our results are valid. In fact, fMRI and ShortPCA both showed expected contralateral cerebral activity around motor regions.
In contrast to Δ Oxy-Hb, our method did not improve estimation using ΔDeoxy-Hb data in Experiment 2 (see details for Supplementary material 4). This could be because the scalp-hemodynamics value is reflected less in ΔDeoxy-Hb and its estimation is difficult owing to the physiological origin of each hemoglobin signal (Haeussinger et al., 2014;Heinzel et al., 2013;Kirilina et al., 2012). Conversely, it could also be considered that cerebral-and scalp-hemodynamics were both poorly detected using Δ Deoxy-Hb because of its low signal-to-noise ratio in our system. If this is the case, and good ΔDeoxy-Hb data is available, our method will be applicable and improvement is expected. Because reduction of scalp-hemodynamics was improved by combining ΔDeoxy-Hb with ΔOxy-Hb (Cui et al., 2010;Yamada et al., 2012), the performance of our method may be further improved by adding ΔDeoxy-Hb information.
The source-detector distance is another methodological concern. We used a source-detector distance of 15 mm to measure the scalphemodynamic artifact because a simulation studies by Okada et al. (1997) reported that the spatial sensitivity profile was confined to the surface layer (i.e., scalp and skull) when the source-detector distance was below 15 mm, and because a 15-mm distance is easier to implement in our current hardware. Although recent studies suggest that signals measured by 15-mm distant Short-channels may contain a small amount of cerebral-hemodynamics, scalp-hemodynamics is likely to be the predominant type of hemodynamics represented in Shortchannel data (Brigadoi and Cooper, 2015;Strangman et al., 2014;Umeyama and Yamada, 2009;Takahashi et al., 2011;Yamada et al., 2009). In a simulation, for example, Yamada et al. (2009) reported that the absorption changes in 15-mm distant channels for gray matter layers were less than 20% of those when the distance was 30 mm. In agreement with these reports, we confirmed by an additional analysis that there were more scalp components in Short-channel data than in Long-channel data in our experimental samples (Supplementary material 9). We evaluated the contribution of scalp-and cerebralhemodynamics for experimental data in Experiments 1A and 2 directly, by VAF. The results revealed that averaged VAFs of the estimated cerebral-and scalp-hemodynamics across samples in the Shortchannels (Experiment 1A) were 0.049 (SD 0.133) and 0.787 (SD 0.122), respectively, whereas those in the Long-channels (Experiment 2) were 0.313 (SD 0.169) and 0.550 (SD 0.163). Thus, scalphemodynamics values were approximately 16 times larger than cerebral-hemodynamics in the Short-channels, but only 1.7 times larger in the Long-channels. Furthermore, the contribution ratio of the 1st PC in Experiment 1A was nearly 0.85 for all samples. These observations strongly support our hypothesis that scalp-hemodynamics is dominant in Short-channels and is distributed globally. Therefore, we expect that multiple 15-mm distant Short-channels are practical enough to accurately extract global scalp-hemodynamics. As our method is applicable to any source-detector distance, future studies can test whether estimation accuracy increases when it is applied to channel distances shorter than 15 mm as has been used in the previous studies (Brigadoi and Cooper, 2015;Gagnon et al., 2012a;Goodwin et al., 2014), provided it can be implemented in the fNIRS measurement hardware.
Another consideration is that because it was designed to remove the global and homogenous artifact, the current method cannot remove local artifacts derived from experimental (motion) and physiological (scalp-hemodynamics) factors. In the motion artifacts, two possible situations may degrade the performance of the proposed method. One is the case in which local artifacts occur in an area in which only Longchannels are placed, and the other is the case in which local artifacts occur in an area in which both Long-and Short-channels are placed. In the first case, local artifacts would only be included in the Longchannels. These local artifacts would not affect the estimation of global scalp-hemodynamics, and would therefore not degrade the performance of other channels. While the global scalp-hemodynamic model would not be able to remove them from those specific Long-channels, they could be removed as a residual when applying the GLM, provided they do not correlate with the task. In the second case, global scalphemodynamics cannot be extracted correctly because of interference by local artifacts included in some of the Short-channels. We simulated both cases with three temporal patterns of motion-like artifacts (boxcar, impulse, and step), which are often observed in fNIRS signals (see Supplementary material 5). In the second case, the contribution ratio of the 1st PC to the Short-channels was markedly lower and estimation error was higher for global scalp-hemodynamics. Therefore, the contribution ratio of the 1st PC could be a good indicator of whether proceeding to the GLM using the estimated global scalp-hemodynamic model is advisable. In fact, when we analyzed data from the patients' right hand movements that had been excluded from analyses because several Long-and Short-channels were contaminated with Step-type local artifacts, we found a much smaller contribution ratio of the 1st PC to the Short-channels (0.539) than was observed in the other data (around 0.85). For these contaminated data, ShortPCA was not more effective than the other methods (see Supplementary material 6). Such local artifacts should be prevented experimentally or removed by preprocessing using other artifact removal methods (e.g., Brigadoi et al., 2014;Cui et al., 2010;Sato et al., 2006).
Even after excluding local motion noise, slight differences among channels remain in scalp-hemodynamics owing to its heterogeneous nature (Gagnon et al., 2012a;Haeussinger et al., 2014;Heinzel et al., 2013;Kirilina et al., 2012). These local components might be eliminated when we allocate the Long-and Short-channels in pairs (Funane et al., 2014;Gagnon et al., 2012aGagnon et al., , 2014Saager et al., 2011;Yamada et al., 2009;Zhang et al., 2009Zhang et al., , 2011. However, Experiment 1-1 showed that more than 85% of the signals measured by the Short-channels contained the global component, with only a small contribution by the local component. Therefore, at least in our experimental setting, using only the 1st PC is effective in estimating cerebral activity. Note that no significant improvement was observed even when we incorporated additional principal components (more than 90% or 100% of contribution ratio; Boden et al., 2007) into the GLM (see Supplementary material 7). As discussed above, large local fluctuations in scalp-hemodynamics (Gagnon et al., 2012a;Kirilina et al., 2012;Yamada et al., 2009) are detectable by a smaller contribution of the 1st PC. In such cases, the proposed method should not be applied.

Advantages of proposed method
By combining four Short-channels with a GLM, our new method improves the estimation accuracy of cerebral activity while keeping the measurement area broad and without reducing the benefits of fNIRS (e.g., few physical constraints and simplicity of use). Hence, this method should be useful, especially in the clinical setting. Consider fNIRS measurement during rehabilitation after stroke (Fujimoto et al., 2014;Obrig, 2014;Kato et al., 2002;Takeda et al., 2007;Vernieri et al., 2006). Unexpected cortical regions, such as those ipsilateral to the moving hand, may be activated during rehabilitation, and detection of these regions requires as broad a measurement area as possible. Furthermore, stroke patients tend to exert maximum effort to move their affected limb, which often causes an increase in scalp-hemodynamics. Using current methods, this results in false positive cerebral activity. Additionally, stress should be minimized as much as possible in clinical practice. Measurement during rehabilitation therefore needs to be from broad cortical areas, with few physical constraints, and with a capacity to remove scalp-hemodynamic artifacts. Previously proposed methods that densely arrange the Short-and Long-channel probes and estimate the local cerebral activity correctly (Funane et al., 2014;Gagnon et al., 2012aGagnon et al., , 2014Saager et al., 2011;Yamada et al., 2009;Zhang et al., 2009Zhang et al., , 2011 do not fulfill these requirements, as probe setting is timeconsuming and the measurement area is limited. Our method meets all the requirements and can be easily implemented in conventional fNIRS systems because the algorithms (PCA and GLM) are simple and the number of Short-channels required is small. Indeed, we confirmed that cerebral activity can be measured without causing a stroke patient significant stress. In addition, we confirmed that our method can also be applied to a verbal fluency task that is commonly used in the clinical setting (Supplementary material 2). Thus, our method is very practical and is expected to be suitable for clinical fNIRS measurement.