2.1.

Subjects and Experimental Setup

The optical data used in this study were obtained from the experimental data previously reported in Franceschini et al.¹¹ The resting fNIRS data have been collected from 17 healthy adults (age range $35\pm 12$ years) lying motionless in a supine position in a quiet dark room. All subjects underwent two to three resting-state sessions of fNIRS measurements, each 300 s. The fNIRS data have been collected with a continuous-wave optical system (CW5 Imager, TechEn Inc., Milford, Massachusetts) comprising 32 laser diodes (average output power of $\sim 1\mathrm{mW}$) at 690 and 830 nm and 16 avalanche photodiodes.¹¹

A whole-head probe has been designed to cover the frontal, parietal, temporal, and occipital lobes of the brain for each subject. Figure 1 shows the schematic of the emitter/detector configuration on the head and the channels defined as the midpoints between pairs of optodes. In total, 50 channels have been used for data collection with a minimum interoptode distance of 3 cm (see Franceschini et al.¹¹ for more details). As shown in Fig. 1, the fNIRS channels were grouped into eight cerebral regions (left/right frontal, temporal, parietal, and occipital regions) for single-subject and group analyses. The optical data were converted to optical density changes as described in Huppert et al.⁹ The modified Beer–Lambert law was used to calculate the concentration change of deoxy-Hb and oxy-Hb based on the differences in the absorption for the two wavelengths. The concentration in total-Hb was calculated by summing the concentrations in deoxy-Hb and oxy-Hb. The average differential path-length factor was assumed to be 6 for both wavelengths.²⁶

Fig. 1

Probe geometry used for data collection.

During fNIRS data recordings, concurrent independent measurements of cardiac activity, respiration, and blood pressure has been made to assess their relationship to spontaneous cerebral hemodynamic fluctuations for each subject. For this purpose, a finger pulse oximeter (Nonin model 8600) has been used to record continuously the arterial saturation and the heart rate from each subject’s left hand. The arterial blood pressure (ABP) has also been measured using a custom-built beat-to-beat noninvasive blood pressure sensor.¹¹ Finally, the respiratory signal has been measured by using a strain gauge respiratory belt (Sleepmate/Newlife Technologies, Resp-EZ). The external systemic physiology signals have been recorded synchronously with the optical data using a National Instrument 6023E card with a sampling rate of 25 kHz. The acquired external data were down sampled to match the sampling rate of the optical measurements (10 Hz).

2.2.

Single-Channel Decomposition of Functional Near-Infrared Spectroscopy Signals

Figure 2 shows the decomposition process comprising three main stages: preprocessing, predecomposition, and single-channel decomposition. In the preprocessing step, the fNIRS signals were first bandpass filtered. Noisy channels were then rejected. Finally, motion artifacts were removed using the spatial ICA. The clean fNIRS signals of all channels were decomposed using the scICA method. We developed a procedure to estimate the optimal number of independent components (sources) for each channel based on the identification of the frequency bands associated with LF oscillations, respiration, and cardiac activity as well as those related to cardio-respiratory interactions. Finally, independent components associated with each frequency band were identified and back-projected onto the measurement space to recover the multichannel fNIRS signals corresponding to that frequency band. These steps will be described in more detail in the following sections. All data processing was done in MATLAB R2011b (The Mathworks, Natick, Massachusetts).

Fig. 2

Schematic of the single-channel decomposition process.

2.2.3.

Estimation of optimal number of sources

For each subject and channel, the estimation procedure was employed in two steps to determine the optimal number of ICs associated with physiological sources.

• Step 1: In the first stage of the estimation procedure, for each subject’s multichannel fNIRS data (deoxy-Hb, oxy-Hb, or total-Hb), we created a frequency-band matrix $M_{\mathrm{FB}}$ with $N_{\mathrm{fb}}$ rows, where $N_{\mathrm{fb}}$ is the number of frequency bands of interest, and two columns including lower- and upper-frequency limits.

  First, an average normalized PSD was computed for each fNIRS channel by decomposing its signal into $K$ independent components using the single-channel decomposition method described in Sec. 2.2.2, where $K$ was the number of the retained principal components (PCs) that accounted for 99% of the variance of the channel’s signal. Then, the PSDs of the $K$ independent components were computed using the Welch spectral estimation method³⁶ with a hamming window of length 1000 samples and an overlap of 20%. The window lengths and overlaps were chosen empirically to achieve high-frequency resolution (0.01 Hz) and good spectral smoothness for the fNIRS data length of 3000 samples (300 s). Finally, the PSDs of the $K$ independent components were normalized to their total power, and then averaged to obtain the average normalized PSD (${\mathrm{PSD}}_{\mathrm{AN}}$) for the channel.

  It can be noted that the normalization procedure was necessary to avoid spectral peaks of different magnitudes masking each other during averaging. The whole procedure was repeated for all $P$ ($\le 50$) channels. Then, a spatially averaged PSD (${\mathrm{PSD}}_{\mathrm{SA}}$) was computed by averaging the ${\mathrm{PSD}}_{\mathrm{AN}}$ of all $P$ channels. By averaging across all channels, we removed spurious spectral peaks of single channels. On the other hand, global spectral peaks with low power were retained through the normalization process. All frequency bands of the ${\mathrm{PSD}}_{\mathrm{SA}}$ were identified by finding the points of inflection using the first and second derivatives. In ${\mathrm{PSD}}_{\mathrm{SA}}$, the spectral peaks (local maxima) were identified by finding the zeros of the first derivative of the spectrum. We then used the zeros of the second derivative of the spectrum to identify the upper and lower limits of each frequency band, i.e., the inflection points on the spectrum, at which the concavity changes from up to down or vice versa. The $N_{\mathrm{fb}}$ spectral peaks and frequency bands were stored in the matrix $M_{\mathrm{FB}}$, which was then used to estimate the number of sources for each channel in the multichannel fNIRS data.

• Step 2: In the second stage of the estimation procedure, for each channel in the fNIRS recording, the optimal number of sources (IC components) was determined between the minimum (${\mathrm{NS}}_{\min }$) and maximum (${\mathrm{NS}}_{\max }$) number of sources. ${\mathrm{NS}}_{\min }$ was set to 3 because in this study, we mainly dealt with three main physiological frequency bands associated with (1) LF fluctuations caused by hemodynamic changes in ABP and local brain hemodynamic activity; (2) respiration (with a spectral peak at $f_R$); and (3) cardiac activity (with a spectral peak at $f_C$). ${\mathrm{NS}}_{\max }$ was determined using the prewhitening procedure by finding the minimum number of PCs ($N_{PC}$) retaining 99% of the fNIRS data variance. Therefore, the optimal number of sources was estimated in the interval [$3N_{PC}$].

  For each fNIRS channel, the estimation procedure started by first decomposing the channel fNIRS signal into three components (${\mathrm{NS}}_{\min }=3$). Then, for each component, a normalized PSD (${\mathrm{PSD}}_{\mathrm{AN}}$) was calculated and analyzed to identify its frequency bands with the method described earlier in this section. The components were grouped based on their frequency content. For this purpose, the number of groups was set to $N_{\mathrm{fb}}$, i.e., the number of frequency bands included in the matrix $M_{\mathrm{FB}}$ of the channel. Therefore, for each frequency band in $M_{\mathrm{FB}}$, all the components whose ${\mathrm{PSD}}_{\mathrm{AN}}$ showed a spectral peak within that frequency band were grouped. Each component was included in multiple groups if its ${\mathrm{PSD}}_{\mathrm{AN}}$ exhibited spectral peaks in all the frequency bands associated with the groups. Then, the grouped ICs were projected back to the original fNIRS signal space and the resulting signal was filtered within the frequency band of the class to remove out-of-band spectral peaks.

  At each run, the residual of the decomposition procedure using $N_S$ sources was computed by subtracting the sum of the reconstructed signals from the original signal. The relative power of the residual was calculated as the ratio of the power of the residual and the total power of the original signal.

  After each run, $N_S$ was increased by an increment of one and the whole estimation procedure was repeated. The optimal number of ICs was finally set to the value with which the relative power of the residual of the decomposition did not change more than 1% by increasing the number of sources (Fig. 3). The main goal of this step was to retain spectral peaks with significant power; therefore, weak spectral peaks, for which no physiological relevance could be found, were not modeled in our study.

Fig. 3

Residual relative power (output of the source number estimation procedure) as a function of number of sources for Subject 1. The curve shows how the residual relative power decreases with increasing number of sources. For this subject, only one fNIRS channel was excluded from the decomposition procedure due to its poor quality. Sources above the upper bound (94) have been rejected by PCA. LB, lower bound; UB, upper bound; and ${\mathrm{NS}}_{\mathrm{opt}}$, optimal number of sources.

2.3.

Single-Subject Spectral Analysis

For each subject, we first determined the association of the spectral peaks in the PSD_SA (see Sec. 2.2.3) of deoxy-Hb, oxy-Hb, and total-Hb with different physiological components. The frequency bands and their spectral peaks with highest power for LF oscillations, respiration ($f_R$), and cardiac pulsations ($f_C$) were identified on the average normalized PSD of the external physiological signals.

The respiratory harmonics were also identified between $f_R$ and $f_C$. We filtered the fNIRS signals above 1.5 Hz; therefore, the cardiac harmonics were not analyzed. To identify frequency peaks associated with nonlinear cardio-respiratory couplings, we used the cross bicoherence method³⁷ to compute the normalized bispectrum of the external respiratory and cardiac signals for each subject. For this purpose, we used the higher-order spectral analysis toolbox in MATLAB.³⁸

For the bispectral analysis, each fNIRS time series was partitioned into 30 nonoverlapping segments considered to be the optimal number of segments required for reliable power estimation.³⁷ The respiration self-coupling and bifrequencies due to interactions between the respiratory and cardiac systems were then identified and labeled on the normalized average PSD of the deoxy-Hb, oxy-Hb, and total-Hb signals.

For each subject, the PSD of the back-projected signals were computed channel by channel and divided into LF fluctuations including very low frequency ($X_{\mathrm{VLF}}$, 0.01 to 0.08 Hz) and LF ($X_{\mathrm{LF}}$, 0.08 to 0.15 Hz) oscillations, respiration ($X_{f_R}$ and its harmonics $X_H$), and cardiac activities ($X_{f_C}$ and cardio-respiratory components $X_{CR}$). Then, for each channel the percentage contribution of these activities were obtained by computing the power ratio of their respective frequency band to the total power of the original signal. The resulting percentage values were averaged over all channels and over those located within different cerebral regions as indicated in Fig. 1. Since each subject underwent two to three resting-state sessions, the results were averaged over the sessions for each subject.

For each subject, we also computed the percentage of residuals resulting from the PCA method ($R_{\mathrm{PCA}}$) and the decomposition procedure ($R_{\mathrm{DP}}$) using the optimal number of sources.

The decomposition was formulated as

In this approach, we assumed that ABP had a global effect across all fNIRS channels.

2.4.

Intersubject Variability

We performed a PCA-based cluster analysis to probe intersubject variability in the percentage contributions of LF fluctuations, respiration, cardiac activity, and other components to the deoxy-Hb, oxy-Hb, and total-Hb signals. For this purpose, we first created a matrix $M$ of size $P\times Q$

3.1.

Spectral Analysis

Figure 4 shows the results of the spectral analysis for Subject 1. The normalized PSDs of the deoxy-Hb, oxy-Hb, and total-Hb signals were computed by averaging the PSDs of all the channels. In Fig. 4(a), three main peaks can be seen at $\sim 0.02\mathrm{Hz}$ (slow activity), 0.13 Hz (respiration), and 1.18 Hz (cardiac pulsations). The magnitudes of these peaks have obscured the other spectral peaks. To extract other spectral peaks, we used the scICA-based decomposition approach described in Sec. 2.2.2. Figure 4(b) shows the normalized PSD resulting from the decomposition of the deoxy-Hb, oxy-Hb, and total-Hb signals. In this subplot, the other spectral peaks are clearly distinguishable from the spectral background. It is of note that in this representation the magnitudes of peaks carry no information because the amplitudes of their corresponding independent components have been normalized to unity during the decomposition process.

Fig. 4

Results of power spectral analysis for Subject 1. (a) Normalized PSD and (b) full spectra averaged over all source–detector pairs for deoxy-Hb, oxy-Hb, and total-Hb signals. To obtain the results in plot (a), the PSD of the deoxy-Hb, oxy-Hb, and total-Hb signals were obtained by averaging the PSDs over all the fNIRS channels. The resulting averaged PSDs were normalized to the maximum power value obtained across deoxy-Hb, oxy-Hb, and total-Hb. To obtain the results in plot (b), each of the fNIRS channels was first decomposed using the single-channel decomposition method (see Sec. 2.2.2). The PSDs of the channel’s independent components were normalized to their total power, and then averaged to obtain the average normalized PSD for the channel. Finally, the averaged PSDs were calculated across channels. In a similar way, plots (c and d) were obtained for the externally recorded ABP, respiration, and cardiac signals for Subject 1.

The same spectral analysis and decomposition procedure were applied to the auxiliary signals recorded synchronously with the fNIRS signals. As shown in Figs. 4(c) and 4(d), the same spectral peaks were found for the auxiliary as those observed for the fNIRS signals. On average across all subjects, the frequency of the systemic ABP varied between 0.05 and 0.15 Hz. The respiration and cardiac activity exhibited a single peak in the frequency ranges [0.13 to 0.36 Hz] and [0.75 to 1.3 Hz], respectively.

As described earlier, we also investigated nonlinear cardio-respiratory couplings using the cross-bicoherence method (Sec. 2.3) for each subject. Figure 5 shows a contour plot of the estimated bicoherence for Subject 1. As shown, the cross-bicoherence method revealed a peak at bifrequency (0.13 to 0.13 Hz) indicating the respiratory self-coupling. The cardio-respiratory coupling has caused a strong bispectral peak at bifrequency (0.13 and 1.18 Hz). All other peaks are related to the nonlinear couplings between the respiratory components (main peak $f_R$, and its harmonics $H_1$ and $H_2$) and the cardiac frequency ($f_C$). The nonlinear interactions appear as $f_C\pm n·f_R$, where $n$ refers to the $n$'th respiratory harmonic. All the peaks labeled at the spectral and bispectral analysis stages were used to compute the percentage contribution of LF, respiratory, and cardiac activities to the fNIRS signals.

Fig. 5

Results of the bispectral analysis between the respiratory and cardiac signals for Subject 1. As shown, the cross-bicoherence analysis reveals phase coupling between the respiratory and cardiac systems. The strong peak at bifrequency (0.13 and 0.13 Hz) indicates the respiratory self-coupling. The peak at bifrequency (0.13 and 1.18 Hz) is related to the cardio-respiratory coupling. All other peaks ($f_2\pm f_1$) belong to the nonlinear couplings between the respiratory components (main peak $f_R$, and its harmonics $H_1$ and $H_2$) and the main peak of the cardiac activity ($f_C$). The color scale is shown normalized to the maximum peak in the image.

3.2.

Contribution of Functional Near-Infrared Spectroscopy Components

On average, the decomposition was performed for all subjects with 37 components and 49 channels after removing noisy channels. Figure 6 shows the average percentage contributions of the LF, respiratory, and cardiac activities and their cross coupling components to the deoxy-Hb, oxy-Hb, and total-Hb signals. For each subject, the averaged percentage contribution of each activity was first computed across all channels for (a) deoxy-Hb, (b) oxy-Hb, and (c) total-Hb. The grand average percentage values were then computed across all subjects. As shown in Figs. 6(a)–6(c), slices 1 and 2, during the resting period two LF components (below 0.15 Hz) were observed. One peaked between 0.01 and 0.08 Hz, associated with slow activities including spontaneous brain hemodynamic activities as well as very LF content of blood pressure variations [shown as VLF in Fig. 6(d)], and the other peaked between 0.08 and 0.15 Hz, mainly related to variations in blood pressure represented as LF in Fig. 6(d). The percentage contributions of VLF were much higher than those of LF in deoxy-Hb, oxy-Hb, and total-Hb.

Fig. 6

Average percentage contribution of LF, respiratory, and cardiac activities to the deoxy-Hb (panel a), oxy-Hb (panel b), and total-Hb (panel c) signals. The pie charts show the contributions of LF activities (slices 1 and 2), respiration (slice 3) and its harmonics (slice 4), cardiac activity (slice 5) and its nonlinear couplings with the respiratory components (slice 6). The percentage of residuals resulted from the PCA method (slice 7) and that from our decomposition method (slice 8) have been also shown in the pie charts (a)–(c). The slices 1 to 6 correspond to peaks 1 to 6 shown in plot (d).

The main spectral peak of the respiration (Fig. 6, slice and peak 3) showed a higher contribution to deoxy-Hb (9.4%) in comparison with the one computed for the spectral peak of the cardiac activity (7.3%; Fig. 6, slice and peak 5 and Table 1). An inverse trend was observed for oxy-Hb and total-Hb. Interestingly, compared to oxy-Hb (7.2%) and total-Hb (5.9%), deoxy-Hb (10.9%) showed higher contributions for the respiratory harmonics (Fig. 6, slice and peak 3) as well as for the nonlinear cardio-respiratory couplings (Fig. 6, slice and peak 6 and Table 1). No significant differences were observed in the residuals of the PCA and decomposition methods between deoxy-Hb, oxy-Hb, and total-Hb [Figs. 6(a)–6(c), slices 7 and 8 and Table 1].

Table 1

Average percentage distribution of sub-components of the slow activity, respiration, and cardiac activities. The percentage values were computed by averaging the values obtained for all subjects. See Fig. 5 for component descripion.

The percentage contributions of the subcomponents of LF, respiration, and cardiac activities have been shown in Table 1. As illustrated, the main contributions belonged to VLF, and the main spectral peaks of respiration and cardiac activities for deoxy-Hb, oxy-Hb, and total-Hb. For slow activities, LF showed higher contributions for oxy-Hb and total-Hb in comparison with deoxy-Hb.

On average across all subjects, in the VLF band, almost 80% of the contributions belonged to the brain autoregulatory activities and the rest to the systemic blood pressure variations. In the LF band, more than 85% of the contributions belonged to the brain autoregulatory activities for deoxy-Hb, oxy-Hb, and total-Hb. Averaged over subjects, the percentage contributions of the VLF components of the brain autoregulatory activity for deoxy-Hb, oxy-Hb, and total-Hb were more than four times the values obtained for the systemic changes in the ABP signal. This ratio was much higher for LF.

Concerning the respiratory subcomponents, the main peak has shown the highest percentage contribution compared to the other subcomponents. However, the first and the second harmonics as well as the nonlinear cardio-respiratory couplings ($f_C\pm f_R$) showed increased contributions to deoxy-Hb in comparison with those found for oxy-Hb and total-Hb. However, the main cardiac peak showed higher contribution to oxy-Hb and total-Hb compared to deoxy-Hb (Table 1).

Figure 7 shows the regional differences in the contributions of LF (sum of VLF and LF), respiration (sum of the main peak and its harmonics), and cardiac activity (sum of the main peak and its nonlinear cardio-respiratory couplings) to the deoxy-Hb, oxy-Hb, and total-Hb signals. To obtain these results, we computed the mean and standard deviation of the contribution values across all channels positioned in the left and right frontal, temporal, parietal, and occipital regions (see Fig. 1) across all subjects. We then statistically compared different cerebral regions using the Mann–Whitney nonparametric test⁴² with $P<0.05$.

Fig. 7

Spatial distribution of contribution of slow activity, respiration, and cardiac activity to the (a) deoxy-Hb, (b) oxy-Hb, and (c) total-Hb signals. The results are represented as the mean and standard deviation of percentage contribution values computed over the fNIRS channels positioned in the left (L) and right (R) frontal (Fr), temporal (Tp), parietal (Pt), and occipital (Op) regions across all subjects.

For deoxy-Hb, statistical comparisons showed significant bilateral decrease in the contribution of the brain slow activity (VLF) moving from the anterior to the posterior parts of the brain. Conversely, the cardiac activity showed higher contributions in the posterior parts of the brain. The right frontal and temporal lobes showed higher percentage contribution for the slow activity with respect to the left side. In all other cases, symmetry was observed between the left and right hemisphere of the brain. In oxy-Hb and total-Hb signals, differences were statistically significant between the occipital regions and frontal areas for LF brain activities and cardiac pulsations. Higher variations were observed in the power of the cardiac activity in the posterior regions across all subjects. For respiration, no significant differences between cerebral regions were observed in deoxy-Hb, oxy-Hb, and total-Hb signals.

3.3.

Intersubject Variability

Figure 8 shows the results of the intersubject clustering. Using the clustering analysis, we identified two distinct physiology types (denoted I and II) included nine and eight subjects, respectively. There were two primary features that distinguished the two physiology types. First, type-I showed a lower contribution of slow activity including VLF and LF components to the deoxy-Hb, oxy-Hb, and total-Hb signals in comparison with type-II. Second, the type-I cluster contained a higher percentage contribution for the respiration (main peak and its harmonics) and cardiac activity (main peak and nonlinear coupling components) in comparison to type-II. Table 2 summarizes all significant differences between physiology type-I and type-II with $P<0.05$.

Fig. 8

Results of the intersubject clustering performed on the percentage contribution of slow activity, respiration, and cardiac activities to the deoxy-Hb, oxy-Hb, and total-Hb signals.

Table 2

Characteristic significant differences between Clusters 1 and 2 (physiology type I and II). VLF, very low frequency; LF, low frequency; fR, respiration, main peak; H, respiration, harmonics; fC, cardiac activity, main peak; fC±n·fR, cardio-respiratory couplings; and NS, not significant.

We further used the unpaired $t$-test with a significance level of $P<0.05$ to determine whether the frequencies of the main peak for slow activity, respiration, and cardiac activities differed significantly between the two subject groups. The subjects exhibiting physiology type-I had significantly lower frequencies for slow activity ($0.012\pm 0.002\mathrm{Hz}$) and respiration ($0.21\pm 0.05\mathrm{Hz}$) in comparison with the second subject group showing physiology type-II with ($0.014\pm 0.004\mathrm{Hz}$) for slow activity and ($0.27\pm 0.05\mathrm{Hz}$) for respiration. No significant differences in cardiac frequencies ($1.07\pm 0.15\mathrm{Hz}$) were found between both physiology types.

For each physiology type, Fig. 9 shows the contribution ratio (oxy-Hb/deoxy-Hb) for slow activity, respiration, and cardiac activity per subject. As shown, for slow activity and respiration, subjects characterized by type-I physiology exhibited slightly higher contribution ratios in comparison with type-II. However, for type-II, the cardiac contribution ratios were greater than those observed for type-I.

Fig. 9

Contribution ratio (oxy-Hb/deoxy-Hb) for slow activity, respiration, and cardiac activities per subject and physiology type.

4.1.

Single-Channel Decomposition

Previous approaches to removing contamination from systemic physiology on fNIRS data have fallen into either spatial or temporal methods. Spatial methods such as PCA^9,43 or spatial ICA methods attempt to remove systemic physiology by projecting out components with the same spatial covariance as systemic noise from the measurements of interest. For example, in the work by Franceschini et al.,¹¹ a resting (baseline) scan was used to construct the spatial components of systemic noise from PCA, which were then projected from the fNIRS scans during the task. However, this approach works well only when the several conditions are met. First, the spatial covariance of the task-related activity must be distinguishable from the background systemic noise. In practice, this requires a large number of fNIRS measurements that cover an area of the head larger then the region activated by the task. Second, the spatial covariance of the physiological noise signal must be stationary. Finally, this approach requires that measurement channels with poor coupling to the head (e.g., high levels of random noise) be removed in preprocessing to accurately estimate the spatial covariance structures. While spatial filtering methods can be applied to analysis of task-related fNIRS signals, these methods cannot be applied to resting-state connectivity studies.

Temporal filtering methods have also been proposed for removing systemic noise. For removal of high-frequency noise such as the main frequency of the cardiac cycle, low-pass linear filtering methods can be applied. However, conventional bandpass filtering will not work properly to remove respiratory, heart-rate variability, or blood pressure-induced fluctuations mainly because the frequency band of interest to be removed is not known a priori and in general are not well separated spectrally from the frequencies of the evoked hemodynamic brain response. Temporal regression methods have been used to remove systemic noise, but require a measure of external physiology with a good quality to remove physiological noise. These methods provide undesirable results when the transfer functions between the sources of activities (here heart and lungs) and the brain are time varying or nonlinear. Diamond et al.⁴⁴ proposed an adaptive regression method using a Kalman filter to address this issue of stationarity of the transfer function from systemic measurements to the brain signal. However, the assumption of a nonstationary linear transfer function from systemic circulation to the brain is an oversimplification of the system. Alternatively, fNIRS measurements with a short separation between the source and detector positions can provide a more regional estimate of systemic physiology since these measurements only probe the superficial skin layers.⁶ These short distances can be used as a local regressor of physiological noise. However, these short-separation signals represent a mixture of different physiological contributions (cardiac, respiratory, and blood pressure) that cannot be separated. When these measurements are used as a linear regressor, an assumption is made that the relative contributions of the physiological signals within the regressor are the same as the contamination in the fNIRS brain measurement. This is only partially true since the deeper measurements contain contamination from both the skin and brain and the brain may have a quite different sampling of arterial and venous contributions from what is captured in the short-distance measurement. When multiple physiological oscillations are contained within a single regressor as is the case in the short-distance correction methods, there is a possibility to introduce artifact frequencies into the measurement of interest, e.g., if ratio of arterial–venous sensitivity differs from the regressor. This is also true in the case of using measurements of systemic physiology from the body. For example, a pulse-oximeter signal from the finger contains a mixture of both cardiac (and harmonics) and respiratory (and harmonics) frequencies. Across the brain, the ratio of these mixed components will vary.

Single-channel ICA is a demixing procedure for decomposing a single time course into its constituent components. Unlike spatial ICA, which decomposes a matrix of data (multiple channels recorded over time) into components, scICA operates on a vector of data (single channel over time) and decomposes the data based on statistically independent components consisting of groups of covarying frequencies from within that one channel. This is an underdetermined decomposition, where a single channel is converted into multiple independent time series. For source analysis in fNIRS signals, we found several advantages for scICA over spatial ICA. First of all, if only a single fNIRS channel signal is recorded or very few channels are available, scICA is able to decompose the signal into its constituent components even with small contribution to the signal.³³ In the case of multichannel recordings, scICA still performs much better than spatial ICA when the activity of interest is highly localized and present only on a few channels. This is because the scICA method utilizes temporal information to isolate underplaying components. Second, the basic assumption to use spatial ICA is that the propagation delays between sources and sensors are negligible. This is not the case for the brain hemodynamic activities with spatially varying delays.⁴⁵ This property affects the identifiability power of the spatial ICA method. However, scICA produces less desirable results under circumstances where the signal generated by underlying sources has components with joint spectral support in the frequency domain.³³

The success of ICA analysis is dependent on several factors. First, the output of spatial and/or temporal ICA depends on the number of components.⁴⁶ If too few components are used, the signal will not be efficiently unmixed. However, if too many components are used, then the underlying signals (physiology) will be split across multiple components. In our approach, the produced components need to be analyzed to identify noise-related components and the relevant ones related to the signal of interest.⁴⁷ To solve these problems, we developed two objective algorithms to efficiently estimate the number of sources as well as to identify the components related to LF oscillations, respiration, and cardiac activity. Our approach for estimating the number of sources searched the optimal number of components between a lower bound, defined as the minimum number of components needed to identify three main frequency bands associated with LF activities, respiration, and cardiac pulsations, and a upper bound, defined as the minimum number of PCs used to reduce data dimensionality. We also used the robustICA method, which has the advantage over the commonly used fastICA algorithm, which is quite sensitive to the initial values.³⁰ In addition, the efficiency of our decomposition approach depended on how well the components related to the respiratory harmonics as well as those related to the nonlinear cardio-respiratory couplings could be identified. To label spectral peaks, we required external physiology. This might be considered as a weakness of our method; however, we only used these reference signals to identify peaks not for the decomposition purpose.

4.2.

Contribution of Functional Near-Infrared Spectroscopy Components

Using scICA, we were able to decompose the deoxy-Hb, oxy-Hb, and total-Hb signals into their underlying components. The benefits of our approach are twofold. First, as we showed, on average across all subjects, 93% of total fNIRS signal power could be decomposed through the decomposition procedure into the constituent components including LF oscillations, respiration and its harmonics, and cardiac pulsations and cardio-respiratory components. The remaining signal power constituted the residuals of our method and PCA. This decomposition helped us to understand the percentage contribution of different physiological sources to deoxy-Hb, oxy-Hb, and total-Hb.

Compared to oxy-Hb and total-Hb, deoxy-Hb showed higher variability across subjects. This is consistent with the results reported by others.⁴⁷ This can be explained by the lower SNR of deoxy-Hb, which made the full decomposition a bit harder. However, during the full spectrum identification process, we found that the powers of fluctuations of different frequency bands in deoxy-Hb were lower than those in oxy-Hb and total-Hb, but their percentage contribution were slightly different from oxy-Hb and total-Hb. Similar observations have been reported by Cui et al.⁴⁸ In our study, as expected, oxy-Hb and total-Hb showed a higher percentage contribution in LF fluctuations and cardiac pulsations in comparison with deoxy-Hb. This is in agreement with the findings reported by other research groups.^5,49

We observed higher contribution (40% to 45%) for very LF fluctuations (0.01 to 0.08 Hz) than that (5% to 10%) found for LF oscillations (0.08 to 0.15 Hz). However, VLF and LF exhibited similar spatial extent.

LF fluctuations have shown a global physiological origin.²⁷ As shown in Fig. 8, we also observed a global effect in VLF and LF in the deoxy-Hb, oxy-Hb, and total-Hb signals. The LF power was more pronounced in the anterior regions; however, symmetry was observed between the left and right hemispheres.

In a recent study, through the analysis of fNIRS data collected in the sensorimotor region, Katura et al.²¹ have shown that 35% of the oxy-Hb signal could be attributed to VLF and LF activities. As shown in Table 1, we found that higher average percentage contribution (50%) for VLF and LF oscillations in deoxy-Hb, oxy-Hb, and total-Hb across all brain areas during the resting state. Duan et al.⁷ suggested that the LF fluctuations in resting-state fNIRS signals may reflect spontaneous brain activity. Since LF fluctuations from background physiology overlap with the LF fluctuations from neuronal activity, we therefore postulate that the increase that we observed in the power of LF fluctuations relates more to the spontaneous background physiology. This may also show that the LF content of fNIRS data decreases in task-induced brain activation.

Using fMRI, a high proportion (75%) of the global LF blood oxygen level dependent noise variance has shown to be attributed to respiratory factors.⁵⁰ In our study, the respiratory main contributions were in the main respiratory peak and its harmonics, which on average constituted 20% of the total power of the deoxy-Hb, oxy-Hb, and total-Hb signals. In the fMRI signal, the respiratory artifact is due to both physiological blood pressure and oxygenation changes similar to those expected to effect the fNIRS signal, but also due to the magnetic susceptibility artifact from expansion of the lungs. This susceptibility artifact only affects the MRI signal and would not show up in the fNIRS and could explain some of the reason why respiratory oscillations account for a larger faction of MRI variance. In addition, fNIRS signals contain contributions from both the skin/scalp and the brain, and this might also account for a discrepancy with the fMRI reports.

4.3.

Variability Between Subjects

On average, we found 20% variability across subjects in different frequency bands. The intersubject variability could be explained by differences between the brain hemodynamic responses of the subjects during the resting state. However, other reasons included the imprecision in the position of the optode between subjects, and removing noisy and artifactual fNIRS channels in the predecomposition stage. In some cases, the latter process reduced the number of channels per brain regions to one. This increased the variability per subject and consequently across subjects.

Our data fell into roughly two types of physiology, with type-I subjects associated with lower relative cardiac and slow component activity, but higher respiratory contributions than type-II. Since this result was for sensors placed all around the head, this finding seems to reflect global physiological patterns and could be explained by differences in the relative contribution of the arterial and venous compartments between subjects in these two groups. Physiologically, this could reflect different levels of baseline cerebral perfusion [e.g., baseline blood flow, volume, vascular volume fractions, or vascular transit time (the ratio of baseline flow and volume; see Ref. 51)] and further investigation is needed to understand these subject-level differences.

4.4.

Application of Single-Channel Independent Component Analysis to Physiological Noise Removal

This work focused on the identification of physiological sources of noise within fNIRS data collected from a probe covering the whole head. A motivation for this study to better understand the nature of this noise is the aim to improve the estimate of brain activity or to characterize neural-related resting-state fluctuations. Using single-channel ICA, we were able to decompose the time course of a single fNIRS measurement into its constituent physiological components. We propose that this decomposition could be used in future work in the context of removing physiology from fNIRS measurements in several ways. First, following the spatial-temporal methods previously described in this context based on conventional PCA or ICA such as those described in Zhang et al.,^9,43 the scICA method could allow more control over removing subcomponents of the signals. In addition, this method could be applied to filter short-separation fNIRS signals wherein the short-separation measurements, which correspond to local scalp measurements, could be further decomposed by scICA and used as multiple regressors. In the short-separation regression methods (e.g., White et al.⁶), the measurement of the scalp signal is used to regress noise from the longer-separation measurements to isolate the brain from the scalp. However, this regression assumes that the relative contributions of the arterial/venous compartments are the same across the probe volume. In other words, a single coefficient is used to model the relationship of the total signal (cardiac, blood pressure, and respiratory) in the short-separation compared to the long-separation measurements. The scICA method applied to the short-separation measurement could provide the means to unmix these physiological components, which could allow multiple regressors to be used in the filtering to account for differences in the relative contributions of these signals across space. Future work is needed to investigate this hypothesis.