1.

Introduction

One way to describe the functional organization of the human brain is to measure resting-state functional connectivity (RSFC) reflecting time synchronized fluctuations in cerebral hemodynamics in the absence of task.^1,2 Spatially distant brain regions showing synchronized hemodynamic activity in these slow spontaneous fluctuations constitute distinct functional patterns of brain activity, also referred to as brain networks.^3,4 The relationship between spontaneous brain activity and task-related function has been observed with paradigms involving different neuroanatomical systems, demonstrating that areas responsible for similar functions also show a correlated activity at rest⁵ and that the response in individual functional areas is correspondingly modulated by the activity of functional resting-state networks.⁶ Moreover, similar patterns of RSFC have been consistently observed across individuals, in both adult⁷ and infant populations.^8,9

RSFC has been mostly investigated with functional magnetic resonance imaging (fMRI) techniques, mainly based not only on the blood oxygenation level dependent (BOLD) contrast but also on measurements of cerebral blood flow with arterial spin labeling,^10,11 or cerebral blood volume with vascular space occupancy imaging.^12,13 The existence of functional brain networks that coexist during RS has also been revealed in direct measurements of neuronal activity with electro- and magnetoelectroencephalography (EEG and MEG),^14,15 suggesting that these signals cannot be attributed only to physiological processes or artifactual components observed in vascular-related functional imaging techniques.

Apart from fMRI and electrophysiological recordings, near-infrared spectroscopy (NIRS) also offers the potential to examine RSFC by measuring concentration changes of oxyhemoglobin (${\mathrm{HbO}}_2$) and deoxyhemoglobin (HbR) in the vasculature of the cortical tissue below measurement channels.^16,17 NIRS has been effectively employed to characterize RSFC in adults¹⁸ and infants,^19–21 and to assess differences in RSFC patterns between experimental groups.²² The most common NIRS RSFC analysis involves evaluating the temporal relationship between time series of the preprocessed data per channel, for example, through relatedness measures, such as the Pearson’s correlation. Individual and/or group inferences are made based on comparing the strengths, directions, and the spatial configurations of the obtained correlation values between groups.

In the current paper, we focus on the specific methodological challenges that could arise during the analyses of RS-NIRS data. The hemodynamic changes elicited by neuronal activity recorded in the NIRS signal contains multiple non-neuronal fluctuations, such as instrumental noise and trends, subject-specific components (e.g., motion-related effects), and hemodynamic fluctuations originated in the cerebral and extracerebral compartments. These fluctuations are related to changes in blood pressure (0.1 Hz), respiration (0.3 to 0.6 Hz), and cardiac pulsation (1 to 2 Hz) and are consistently present at specific frequencies and across measurement channels.²³ The existence of these non-neuronal components introduces spurious common variance across time series, and therefore could falsely increase FC between signals. Furthermore, changes in ${\mathrm{HbO}}_2$ and HbR signals are also affected by the filtering effect of the actual hemodynamic response function (HRF) that acts as a low pass filter with cutoff frequency approximately around 0.1 to 0.2 Hz. Consequently, the NIRS signal has a colored spectrum (i.e., it is not a white process with equal energy across all its frequencies), and it exhibits nonzero temporal autocorrelation. This fact is relevant when RSFC is examined based on pairwise correlation measures. As first noted by Granger and Newbold,²⁴ the correlation between two random signals will artificially increase if these signals exhibit nonzero autocorrelation. Therefore, the intrinsic autocorrelation in the NIRS data may artificially inflate correlation values, increase the false-positive rate under the null hypothesis of no correlation between channels, and potentially compromise the validity of data analyses outcomes.^25–28

The impact of autocorrelation in the validity of statistical estimation has been widely discussed in the analysis of task-related activity in fMRI data^29–31 but only recently in the NIRS literature.^25,27,32,33 The most common way to account for signal autocorrelations (also referred to as serial correlations) is to prewhiten the signal so that the residuals of a linear regression model, which describe the hypothesized task-related activity, become uncorrelated (i.e., white).^25,29–32 The use of prewhitening in the analysis of RS data is less common and only few studies have suggested prewhitening of the signal itself, either in fMRI data,^34,35 or NIRS data.²⁶ By prewhitening RS data, it is expected that the signals become white, exhibiting less temporal autocorrelation, thus reducing false-positive rates. Yet, this approach also presents some caveats. We might assume that the colored frequency spectrum of the RS-NIRS signal is partially originated by physiological noise or induced by the blurring of the hemodynamics, which act as a low pass filter. However, spontaneous neural oscillations as measured with electrophysiological techniques also exhibit a colored frequency profile, typically characterized as a $1/f$ spectrum.^36–38 In this scenario, convolving a neural signal with $1/f$ frequency profile with the HRF and adding physiological noise, as is the case in RS-NIRS signals, would still leave a signal that has a $1/f$ profile in the limited band of the HRF. Therefore, it might not be a reasonable null hypothesis to expect RS-NIRS signals to be white, as this would imply that the information related to the HRF has been partially or completely removed, which would complicate the neurophysiological interpretation of the signal that is left after applying prewhitening.

Because prewhitening has been only recently proposed as an important step in the RS-NIRS preprocessing pipeline,²⁶ a careful examination of its effect on the physiological properties of the RS-NIRS signal for FC analyses is necessary. Therefore, the purpose of the current work is to evaluate the effect of prewhitening during the analyses of real RS-NIRS data. For this purpose, we replicated two previous NIRS studies that measured RS activity in infants, by following similar data acquisition, data preprocessing, and analyses procedures.^19,21 Also, the results obtained following a standard preprocessing pipeline (as done in the previous studies) are compared with the results after incorporating two common prewhitening approaches in the preprocessing pipeline. The two prewhitening approaches include one based on a nonparametric approach, i.e., assuming no model for the autocorrelation as it is used for fMRI data analyses³¹ and one based on modeling the signal as a stochastic autoregressive (AR) process.²⁵ By comparing these two algorithms, we also assess whether different prewhitening procedures produce different results. The analyses described above were performed on a data set collected from 4-month-old sleeping infants, obtained in a similar testing procedure to those of Homae et al.¹⁹ and to the recent study by Watanabe et al.²¹ All evaluations were performed for both ${\mathrm{HbO}}_2$ and HbR.

2.

Methods

2.1.

Data Acquisition

Data from 24 healthy infants were included in this study (age: $124.6\pm 3.76$ days, 12 female). Parents were informed about the procedure and then signed a written informed consent before the experiment. This study was approved by the local ethical committee. Spontaneous hemodynamic activity was recorded using NIRS with a NIRScout system (NIRx Medical Technologies, California) at wavelengths 760 and 850 nm. About 16 light emitters and 24 detectors were positioned over frontal, temporal, parietal, and occipital regions of both hemispheres according to the international 10–20 system, using the nasion, inion, and preauricular points as external head landmarks (Fig. 1). Every pair of adjacent light emitters and detectors formed a single measurement channel, which allowed for simultaneous recordings of 52 channels for each hemoglobin oxygenation state at a sampling frequency of 8.93 Hz. A diagram showing the distance between adjacent source–detector pairs used as measurement channels is presented in Fig. 1. The average distance between adjacent source–detector pairs in this study was 2.5 cm. Occipital channels were discarded for all participants for being particularly prone to motion artifacts. Thus, only data from the remaining 14 sources and 19 detectors (46 channels) was considered in the analyses. During data acquisition, participants were in natural sleep on their caregivers’ lap in a sound attenuated room, without another source of illumination than the screen of the recording computer. Parents were asked to remain silent and to minimize movements over the duration of the recordings in order to avoid involuntary cap or optode displacement. Recordings lasted between 11 and 21 min.

Fig. 1

Diagram of fNIRS optode (sources = red rhombus, detectors = gray circles) and channel (numbers) localization in the cap employed in the current experiment. The right part of the figure shows a diagram with the distances between sources and detectors.

2.3.

Data Preprocessing

First, light intensity data (raw data) were converted into optical density changes (hmrIntensity2OD in Homer2). Noisy segments typically occurring at the beginning and at the end of each dataset, corresponding to awake activity of the infants (before the infant fell asleep and after the infant woke up), were visually identified and discarded. Although most participants displayed good data quality, some datasets contained motion induced spikes and signal drifts that were corrected using the wavelet-based despiking method described in Patel et al.,⁴⁰ which was adapted for the NIRS data. This method effectively reduces motion artifacts in the data while keeping time periods free of artifacts unaffected (see Fig. 2).

Fig. 2

Time-series of a representative example before and after wavelet despike, and removed noise at this step. This wavelet-based method⁴⁰ is designed to detect and remove transient events caused by outliers in the signal and low-frequency trends. In the top part of the right column of figure, the adjacency matrix of the clean data after conversion to ${\mathrm{HbO}}_2$ and HbR and bandpass filtering is shown. The plot in the bottom part of the right column shows the power spectrum of the clean data time-series. This figure was created for every participant as part of the data quality assessment routine during preprocessing.

Subsequently, optical density data were converted into ${\mathrm{HbO}}_2$ and HbR concentration changes (hmrOD2Conc in Homer2), considered as differential path length factors 5.3 (760 nm) and 4.2 (850 nm).⁴¹ The signal was bandpass filtered between 0.005 and 1 Hz to remove the contribution of very slow frequency fluctuations, as well as high frequency components (e.g., cardiac pulsation).¹⁹ After this step, all datasets were limited to $\sim 571\mathrm{s}$ (i.e., 5100 samples) to ensure a homogenous contribution across participants to the statistical analyses. This step was performed by visually inspecting the data in order to select the segment displaying the best data quality. The first 150 samples and the last 5 samples were also removed to avoid the beginning and end effects originated by the prewhitening methods. Thus, the final duration of the time series was $\sim 554\mathrm{s}$ (i.e., 4945 samples).

2.7.

Simulations

Two sets of simulations were performed to further investigate the effect of prewhitening in signals exhibiting temporal autocorrelation and to facilitate the interpretation of the results with real RS-NIRS data.

2.7.1.

False-positive rate in randomly generated time series

To investigate the effect of prewhitening on the Pearson correlation coefficients, two sets of randomly generated time-series of length of 5000 sample points were simulated following a similar procedure, as described in Santosa et al.²⁶ One set of simulated time-series was generated by filtering two random normally distributed signals with a $1/f$ shape filter, in order to make these signals match the frequency spectrum of spontaneous neural oscillations. The second set of simulations was generated with a white frequency spectrum. In addition, two different conditions of correlation were assessed for each set: with no correlation (i.e., the distribution of the correlation coefficients is expected to be centered in zero) and with an induced correlation of $r=0.5$, which was generated by multiplying the time-series with an upper triangular matrix obtained by the Cholesky decomposition of the desired correlation matrix (i.e., [1 0.5; 0.5 1]).

The following processing steps were independently applied in each of the four simulated datasets. First, signals were convolved with an HRF generated using a gamma function $\mathrm{\Gamma }(t;n,\lambda )$ with parameters $n=4$ and $\lambda =2$. Then, each signal was added terms simulating physiological noise (i.e., cardiac pulse and respiration) that were generated as narrow-band filtered noise (0.25 to 0.35 Hz for respiration, 1.95 to 2.05 Hz and 3.95 to 4.05 Hz for cardiac pulse), as well as a white noise component that simulated hardware-related random noise. The amplitude of these components relative to the HRF signals was 1/2 for the cardiac pulse, 1/20 for the respiration, and 1/10,000 for the hardware-related random noise, which yielded signals with similar power spectra to our experimental data (see Fig. 3). Following the same preprocessing pipeline as for the experimental data, signals were bandpass filtered between 0.005 and 1 Hz. Finally, prewhitening was applied independently to each signal by the two procedures described above (i.e., nonparametric approach and AR fitting approach). In each repetition ($n=3000$), Pearson’s correlation for each pair of simulated signals was computed after every preprocessing step to evaluate the effect of each of these steps in the distribution of correlation coefficients.

Fig. 3

Example of a simulated signal and its power spectrum at each processing step, and power spectrum of ${\mathrm{HbO}}_2$ for a representative example of our experimental data.

2.7.2.

Effect of prewhitening on simulated RSFC patterns

To further investigate the effect of prewhitening on the analysis of RSFC, 20 datasets with “RSFC patterns” similar to those observed in real experimental data were simulated following the same procedure described above. Each dataset included 40 time-series of 5000 sample points, thus each of these datasets simulated RS data of a participant acquired in a similar testing procedure as our experimental data. The original time-series were random uncorrelated signals with a $1/f$ frequency spectrum, which could be regarded as signals reflecting spontaneous neural processes. Different levels of correlation between these simulated “channels” were induced, producing adjacency matrices showing specific and known correlation patterns. The same correlation structure was imposed in the 20 simulated “participants.” The correlation structure was induced to approximately match the shape of the adjacency matrices observed in our experimental RS-NIRS data. A correlation of $r=0.8$ was induced between neighboring channels, whereas a correlation of $r=0.6$ was induced between homotopic channels. Correlations of $r=0.1$ or $r=0.2$ were simulated between the remaining pairs of channels. Similarly, these signals were convolved with the HRF, corrupted with physiological noise to create simulated NIRS datasets and then bandpass filtered. The two prewhitening procedures were then applied to these datasets. Changes in the spatial patterns of correlation were investigated by looking at the group level adjacency matrices of the simulated NIRS data, after bandpass filtering and after prewhitening. The same hierarchical clustering approach as described above was also performed, to investigate whether the induced spatial clustering structure was preserved after the different processing steps.

3.

Results

3.1.

Functional Connectivity Analyses

Figure 4 illustrates the effect of the different preprocessing steps and prewhitening methods on the ${\mathrm{HbO}}_2$ and HbR signals [Fig. 4(a)] and the corresponding initial 300 autocorrelation coefficients [Fig. 4(b)] for a single channel of an exemplar subject. Without bandpass filtering, the HbR signal shows a high autocorrelation value even after 300 time lags ($\sim 30\mathrm{s}$), whereas the autocorrelation coefficients of the ${\mathrm{HbO}}_2$ signal become approximately zero after $\sim 200$ time lags ($\sim 20\mathrm{s}$). After band-pass filtering (0.005 to 1 Hz), the ${\mathrm{HbO}}_2$ and HbR signals show similar autocorrelation that decreases to zero after $\sim 200$ time lags. As shown in the bottom plots of Fig. 4, both nonparametric and AR prewhitening procedures reduce the autocorrelation of the data, even though the autocorrelation coefficients shrink to values closer to zero at shorter lags using the AR approach than the nonparametric approach.

Fig. 4

(a) The times-series of a single channel for ${\mathrm{HbO}}_2$ (red) and HbR (blue) at different steps of the preprocessing. (b) The first 300 autocorrelation coefficients.

Figure 5 shows the FC matrices for the three preprocessing approaches for an individual subject (first row) and at the group level (second row). It is generally assumed that ${\mathrm{HbO}}_2$ and HbR signals in NIRS data should exhibit negative correlation.^21,43 This assumption is fulfilled both at the individual and at the group level when data has only been preprocessed with band-pass filtering. Moreover, homologous regions of both hemispheres show high correlation values in ${\mathrm{HbO}}_2$ and HbR, forming clusters and a clear spatial distribution around the main diagonals. In contrast, incorporating any of the two prewhitening procedures as part of the preprocessing makes the negative correlation between ${\mathrm{HbO}}_2$ and HbR disappear, and the correlation between homologous regions is only evident in ${\mathrm{HbO}}_2$, with a different spatial distribution in each of the methods.

Fig. 5

FC matrices for an individual subject (first row) and at the group level (second row) for the three types of preprocessing. In each plot, the FC matrix for ${\mathrm{HbO}}_2$ and HbR is shown in the top-left part and bottom-right part, respectively (RH: channels in right hemisphere, LH: channels in left hemisphere). The FC matrix representing the correlation between ${\mathrm{HbO}}_2$ and HbR is shown in the top-right. Note that the plots are symmetric with respect to the main diagonal.

We further investigated changes in correlation values due to the effect of prewhitening across ${\mathrm{HbO}}_2$, HbR, and between ${\mathrm{HbO}}_2$ and HbR (Fig. 6). For that, correlation coefficients were transformed to $z$ scores by Fischer’s z transformation. For each condition (i.e., ${\mathrm{HbO}}_2$, HbR, and ${\mathrm{HbO}}_2\text{-}\mathrm{HbR}$), the individual $z$ scores were evaluated as random effects and ANOVA was performed with preprocessing method as a factor (3 groups of 24 infants in total). Posthoc pairwise multiple comparison tests were performed following the Tukey’s honestly significant difference criterion (Table 1). For ${\mathrm{HbO}}_2$, we observed a significant difference between the three methods ($F_{[\mathrm{2,69}]}=8.62$, $p=0.0004$). Multiple comparison test showed that correlation values in ${\mathrm{HbO}}_2$ for the AR method were larger than those of the nonparametric approach ($p=0.0002$). For HbR, we also found a significant difference between the methods ($F_{[\mathrm{2,69}]}=146.13$, $p<0.00001$). Posthoc tests revealed that $z$ scores were larger for the bandpass filter method than those of the nonparametric approach ($p<0.00001$) and the AR approach ($p<0.00001$). For the ${\mathrm{HbO}}_2\text{-}\mathrm{HbR}$ condition, a significant difference between the methods was also observed ($F_{[\mathrm{2,69}]}=91.9$, $p<0.0001$). Pairwise posthoc tests demonstrated that $z$ scores in the bandpass filter condition were smaller (i.e., larger negative correlation values) than those of the nonparametric approach ($p<0.00001$) and those of the AR approach ($p<0.00001$).

Fig. 6

Histograms of the group level correlation $r$ values for the different preprocessing methods. Each column shows the comparison between preprocessing methods for the different conditions (i.e., ${\mathrm{HbO}}_2$, HbR, and ${\mathrm{HbO}}_2-\mathrm{HbR}$). Each infant contributed with $46\times 45/2=1035$ correlation $r$ values (upper triangular connectivity matrix) to each histogram.

Table 1

Mean z score differences between preprocessing methods and 95% confidence interval.

Method	Mean difference z	95% Confidence interval
${\mathrm{HbO}}_2$
Bandpass filter—nonparametric	0.0662	$-0.0073$ to 0.1397
Bandpass filter—AR	$-0.0612$	$-0.1347$ to 0.0123
Nonparametric—AR	$-0.1274$	$-0.2009$ to $-0.0539$
HbR
Bandpass filter—nonparametric	0.2904	0.2403 to 0.3405
Bandpass filter—AR	0.3254	0.2754 to 0.3755
Nonparametric—AR	0.0350	$-0.0150$ to 0.0851
${\mathrm{HbO}}_2\text{-}\mathrm{HbR}$
Bandpass filter—nonparametric	$-0.2464$	$-0.2927$ to $-0.2000$
Bandpass filter—AR	$-0.2015$	$-0.2479$ to $-0.1552$
Nonparametric—AR	0.0448	$-0.0015$ to 0.0912

3.2.

Hierarchical Clustering Analysis

The results of the hierarchical clustering analysis of the standard preprocessing were consistent with the results previously reported in Homae et al.¹⁹ A large degree of reproducibility is observed despite differences in data acquisition between studies (46 channels as opposed to 94 channels) and recording duration for each infant (9 min instead of 3 min). Also, we assessed 4-month-old infants, which correspond to an intermediate age between the infants at 3 and 6 months of age tested in the previous study.¹⁹

As it can be seen in Fig. 7 for ${\mathrm{HbO}}_2$, if only band-pass filtering is applied, clusters are formed between homologous regions of both hemispheres, which split into individual frontal, temporal, and parietal regions on each hemisphere at a lower threshold (i.e., larger degree of similarity). With the nonparametric and AR prewhitening, the clustering between homologous regions is preserved, but each method shows a different spatial pattern. With the nonparametric approach, the most posterior parietal channels in both hemispheres are clustered along with the most frontal channels. The AR prewhitening approach shows a similar distribution to the nonparametric approach, except the anterior–posterior cluster is formed between channels located in the frontotemporal region and the posterior parietal region. For the lower threshold, AR prewhitening results into a larger number of interhemispheric clusters than the nonparametric approach, being the clustering in the later condition more similar to that obtained with band-pass filtering.

Fig. 7

Hierarchical clustering analysis for ${\mathrm{HbO}}_2$ and HbR data. Top two rows show results for the higher threshold corresponding to three clusters. The two rows at the bottom of the figure show results for the lower threshold corresponding to six clusters.

As for HbR data, the spatial distribution of clusters with only band-pass filtering at the higher clustering threshold (i.e., few clusters) is comparable to that obtained in ${\mathrm{HbO}}_2$ and to results presented in a similar study by Homae et al.⁴⁴ The results at the lower clustering threshold (i.e., more clusters) and with the two prewhitening procedures show a less bilateral spatial distribution, which is restricted to a few channels. A frontal cluster involving both hemispheres is formed in the three approaches, regardless of the thresholds, but the main temporal/parietal clusters appear separated across hemispheres.

3.3.

Phase Difference Analysis

Figure 8 shows the effect of prewhitening in hPod values across channels for a representative subject [Fig. 8(a)] as well as at the group level [Figs. 8(b) and 8(c)]. After band-pass filtering, we observed a pattern of hPod values that follows previous observations for a group of participants with a similar age range as our data.²¹ In both plots of the standard and polar histograms with only band-pass filtering, it can be seen that most hPod values are close to $\pi$ rad (i.e., 180 deg), which corresponds to an antiphasic pattern between both ${\mathrm{HbO}}_2$ and HbR signals. However, prewhitening disrupts this pattern. The overall distribution of hPod values after prewhitening becomes more widespread, although the distribution for prewhitening with the AR approach shows a mode closer to $\pi$ (i.e., antiphase) than the nonparametric approach. By using a Watson–Williams test, we observed a statistically significant difference in average hPod values between the methods ($F_{[\mathrm{2,69}]}=63.4342$, $p<0.0001$). Posthoc pairwise comparison tests between the methods revealed that mean hPod values in the bandpass filter method were larger than those of the nonparametric approach ($p<0.0001$) and the AR approach ($p<0.0001$). Significant differences were also observed between the two prewhitening methods, the AR approach showing larger hPod values than the nonparametric approach ($p=0.0004$; Tables 2 and 3). Note that values range between 0 and $\pm \pi$, so higher values correspond to values that are closer to $\pm \pi$.

Fig. 8

Polar histogram showing phase difference (hPod) values (a) for an individual subject and (b and c) at the group level for the three types of preprocessing. Histograms in the panel (c) show absolute hPod values.

Table 2

Group average hPod values and standard deviation for each preprocessing method.

Method	hPod	Standard deviation
Bandpass filter	$-2.8668$	0.36
Nonparametric	$-1.8342$	0.91
AR	$-2.2697$	0.88

Table 3

Pairwise comparisons on mean hPod difference values between preprocessing methods.

Method	F[1,46]	P value
Bandpass filter—nonparametric	144.840	$<0.0001$
Bandpass filter—AR	64.6894	$<0.0001$
Nonparametric—AR	14.1282	0.0004

3.4.

Simulations

Figure 9 depicts the histograms of Pearson correlation coefficients after each processing step for each simulated dataset, where the signals with $1/f$ spectrum and white spectrum are shown on the top and the bottom, respectively. For the two sets of simulated data without correlation ($r=0$), the distribution of $r$ coefficients becomes wider after convolution with the HRF and bandpass filtering, suggesting that both steps artificially increased correlation values and consequently increased false-positive rates. After prewhitening, correlation coefficients become centered in zero again, following a similar distribution as the original signal with the nonparametric approach, whereas the distribution is narrower with the AR approach. For induced correlation $r=0.5$, the distribution of correlation coefficients is wider for the signals with a $1/f$ power spectrum than for the white signals due to their autocorrelation. Similar to the scenarios with $r=0$, convolving with the HRF widens the distribution of correlations, and bandpass filtering reduces the confounding effects due to physiological noise. Decisively, both prewhitening approaches substantially reduced the correlation between the time-series towards zero. The nonparametric prewhitening resulted in distributions centered in zero for datasets with $1/f$ spectrum in the two correlation conditions. Similar results were obtained for the dataset with white spectrum and $r=0$. The dataset with a white spectrum and $r=0.5$ showed a distribution slightly shifted toward positive correlation values. The AR approach showed a similar distribution, with a narrow shape and centered above zero, across datasets and conditions.

Fig. 9

Distribution of the measured Pearson correlation coefficients for two sets of simulated time-series under different levels of correlation. The distribution of correlation coefficients is plotted for each processing step. (a) Uncorrelated signals with $1/f$ power spectrum. (b) Signals with $1/f$ power spectrum and an induced correlation of $r=0.5$. (c) Uncorrelated signals with a white power spectrum. (d) Signals with a white power spectrum and an induced correlation of $r=0.5$.

Figure 10 shows the group level adjacency matrices for the four conditions assessed in the second set of simulations. Adjacency matrices for the simulated neural data and for the NIRS data show the expected induced correlation patterns. By contrast, the two prewhitening procedures reduce the correlation between channels toward values close to zero, in agreement with the results shown in Fig. 9. Prewhitening also alters the spatial configuration of the adjacency matrices, although part of the induced correlation structure is preserved (note differences on the scale of the plots). Figure 11 illustrates that the hierarchical clustering analysis for the simulated neural data and for the NIRS data shows the expected spatial clustering configuration with the “anterior” and “posterior” channels clustering together and splitting by “hemisphere” when the threshold is set to display larger degrees of similarity (i.e., dark and light for blue and red colors). After prewhitening, the “anterior” and “posterior” clusters are preserved at the higher threshold, but the clustering configuration is disrupted when clusters between more similar channels are displayed. With nonparametric prewhitening, the spatial configuration changes and clusters that do not belong to the original structure are formed. The AR prewhitening approach also results in diffferent clusters from those of the original configuration.

Fig. 10

Group level adjacency matrices from simulated data showing specific RSFC patterns resembling real data. (a and b) The results for the simulated neural and the simulated NIRS data. (c and d) The adjacency matrices for the two prewhitening procedures. Note that the scale in each row is different.

Fig. 11

Results of the hierarchical clustering analysis corresponding to simulated “participants.” (a and b) The results for the simulated neural and NIRS data. (c and d) The results after applying each prewhitening procedure. Each color represents channels that are grouped together forming clusters based on similarity. For the blue and red clusters, dark and light colors indicate that these channels belong to the same cluster if a higher threshold is considered, but split into two clusters at lower thresholds (i.e., larger degree of similarity).

4.

Discussion

It has been previously described that the validity of the standard NIRS RSFC data analysis approach is compromised by the presence of physiological and artifactual fluctuations that increase the level of autocorrelation in the data as well as introduce shared variance across different channels, leading to uncontrolled false-positive rates.²⁶ One of the main causes of this issue is the presence of autocorrelation in the RS-NIRS data due to the slow nature of hemodynamic fluctuations in the extracerebral and cerebral compartments.²⁷ To overcome this problem, it has been proposed that autocorrelation should be removed from the data before performing RSFC analyses via prewhitening procedures.^26,34,35

Here, using real RS-NIRS data obtained from infants, we assessed the effects of different prewhitening methods, in comparison to a standard analysis method (e.g., one excluding prewhitening). The rationale presented in Santosa et al.²⁶ to recommend the use of prewhitening in RS-NIRS data, which has also been proposed for RS-fMRI studies, ^28,34 is based on the fact that temporal autocorrelation inflates relatedness measures, such as Pearson’s correlation. In their work, Santosa et al.²⁶ convolved two random uncorrelated signals with a canonical HRF and added physiological noise (i.e., cardiac pulse and respiration) to make these signals show the properties of real RS-NIRS signals. They demonstrated that this step make signals self-correlated, increasing the correlation between them and producing false-positive rates around 80%. In their work, prewhitening indeed removed temporal autocorrelations, and false-positive rates returned to the expected true value of 5%. Nevertheless, it is also true that real hemodynamic signals are not expected to show a white frequency spectrum. Therefore, prewhitening the NIRS signal and consequently altering the expected power law behavior of RS-NIRS data in the frequency band of the HRF may compromise the interpretation of RSFC results (and its link with previous literature assessing RSFC), as the prewhitened signals will no longer hold the specific properties of spontaneous hemodynamic fluctuations as measured with fMRI and fNIRS (i.e., BOLD or ${\mathrm{HbO}}_2$ and HbR).

In the current work, we specifically evaluated the effect of the following three preprocessing methods in experimental RS-NIRS data obtained from 4-month-old infants: (1) the “standard” band-pass filtering approach, (2) prewhitening based on a nonparametric algorithm, and (3) prewhitening based on an AR fitting approach ²⁵. We also compared our results to previous studies assessing RS activity in infants using NIRS.^19,21 We note that our two prewhitening procedures succeeded in reducing the autocorrelation of the signal. AR prewhitening was implemented following Barker et al.³² and Santosa et al.²⁶ The range of the optimal AR model in our study, which varied from $P=60$ up to $P=110$ at a sampling rate of 8.93 Hz (i.e., a range from 6.7 to 12.3 s) is similar to those observed in Santosa et al.,²⁶ which reported model orders up to $P=40$ at a sampling rate of 4 Hz (i.e., also $\sim 10\mathrm{s}$). Apart from the sampling rate, the order of the optimal AR model also depends on signal quality and varies across measurement channels.

Our study yielded three relevant findings. First, we were able to replicate the results of two previous infant RS-NIRS studies ^19,21 when using the traditional preprocessing pipeline including band-pass filtering but not prewhitening. Second, results after data prewhitening were different from those presented in the original studies. Third, the two different prewhitening methods produced different results. Considering the hierarchical clustering analysis of the ${\mathrm{HbO}}_2$ signal for both prewhitening methods (Fig. 7), a cluster comprising channels located in frontal and posterior parietal regions emerged that is neither observed in the results without prewhitening nor in the original results.¹⁹ Most networks/clusters reported in the infant RS literature^8,19,44–46 are formed between bilateral regions (i.e., frontal, temporal, or visual), but there are also reports of a developing frontoparietal network,⁴⁷ which has been consistently reported in studies with adults. This frontoparietal network disappears when the threshold is defined to reveal clusters showing larger similarity (i.e., a larger number of clusters). Spatial clustering for HbR occurs between channels proximate to each other and less interhemispheric clustering is observed, with similar results across the two prewhitening methods. Our results differ from the results shown in Homae et al.⁴⁴ for HbR, which were characterized by interhemispheric clustering with almost identical distribution as ${\mathrm{HbO}}_2$. To our knowledge, Homae et al.⁴⁴ is the only study that reported HbR with RS-NIRS infant data to date.

Perhaps, one of the most relevant observations of our analyses is that the negative correlation and antiphase state between ${\mathrm{HbO}}_2$ and HbR signals was no longer visible after prewhitening (Figs. 5, 6, and 8), although there is a slight tendency toward an antiphase relationship with the AR approach. These results contradict previous observations of the intrinsic physiological relationship between ${\mathrm{HbO}}_2$ and HbR in the NIRS signal, ^16,21,43,48 which we believe is an important discussion point. The signal autocorrelation originates from all its sources, including not only artifactual and non-neurobiological physiological confounds but also “brain-related” hemodynamic activity, as the HRF acts as a low-pass filter.²⁶ Therefore, removing autocorrelation might imply that the HRF of interest is also partially or completely removed, and this might cause the expected negative correlation between ${\mathrm{HbO}}_2$ and HbR signals to disappear. This opens a question on the neurobiological significance of the remaining signal (i.e., the link between neural processes and the prewhitened signal) and raises concerns about the effects of using prewhitening in RS-NIRS studies. We could assume that signals after prewhitening (i.e., innovations) do not represent ${\mathrm{HbO}}_2$ and HbR fluctuations, but they reflect the activity-inducing signals underlying them. From that perspective, prewhitening could be understood as a hemodynamic deconvolution procedure, as it has been proposed for fMRI data analysis.^49–51 However, this is not the primary goal of prewhitening as it is understood in this work (i.e., removing the signal autocorrelation). In our opinion, and on the light of these results, further investigation on the neurobiological relevance of the prewhitened signal in RS-NIRS data is still required before directly adopting this procedure.

For both prewhitening approaches, we observed a significant reduction in the correlation values obtained in the adjacency matrices of HbR signals but not in those of ${\mathrm{HbO}}_2$ signals (Figs. 5 and 6). We note that motion induced artifacts, although sparse, might still be present in our experimental RS-NIRS data, and potentially modulate the correlation between channels. However, we visually inspected prewhitened data and observed that time-points presumably showing motion related artifacts co-occurred across ${\mathrm{HbO}}_2$ and HbR. Therefore, the differences on distribution of correlation values after prewhitening observed across ${\mathrm{HbO}}_2$ and HbR cannot be attributed to data showing different motion-related effects across components. Alternatively, we infer that these results could be explained by the fact that ${\mathrm{HbO}}_2$ is more sensitive to physiological fluctuations than HbR, and thus be more prone to nonstationary low-frequency physiological fluctuations related to modulations in cardiac rate,⁵² breathing,⁵³ and systemic blood pressure.^54,55 These processes are difficult to be explained by means of stationary models of correlations but could modulate time–frequency dynamics of hemodynamic fluctuations during resting-state, and account for a substantial part of the correlation between the channels in our experimental ${\mathrm{HbO}}_2$ data after prewhitening.

Finally, our simulations demonstrated that convolving two uncorrelated time series with the HRF artificially increases the correlation between them, resulting in larger false-positive rates. These results were already noted by Santosa et al.²⁶ in similar simulations, and they are also in line with those reported in Bright et al.²⁸ for simulations with fMRI-like data. Figure 2 in Bright et al.²⁸ illustrated that the false-positive rate increases due to the effect of bandpass filtering, specifically that narrower bandpass filters and faster sampling rates show the highest false-positive rates. Thus, the effect of bandpass filtering the data can be considered equivalent to convolution with the HRF. In all the conditions (i.e., $r=0$ and $r=0.5$), AR prewhitening removed correlation between the signals and centered the distribution of the measured correlation $r$ values in zero. Nonparametric prewhitening also removed the correlation and centered the correlation $r$ values in zero, except for the signals with a white spectrum and $r=0.5$, where the distribution of correlation coefficients was slightly shifted toward positive values.

In the second set of simulations, where changes in spatial correlation patterns resembling RSFC were investigated, we observed similar correlation patterns between the “original data” (i.e., simulated neural activity where the correlation structure was induced) and the simulated “NIRS data.” After prewhitening, and with the two approaches, correlation between channels was reduced to levels close to zero, and the shape of the adjacency matrices disrupted, although part of the original correlation structure is still preserved. It is important to note that in our simulations, all the “participants” showed the same original correlation structure, which we consider a very unlikely scenario in the case of real RS-NIRS data. We expect that, when experimental data are considered, RSFC results of the measured RS-NIRS data will considerably differ from results after applying prewhitening.

5.

Conclusion

This work investigated the implications of including prewhitening as part of the RS-NIRS data preprocessing pipeline. More specifically, the effect of prewhitening was assessed by replicating the results of two previous RS-NIRS studies and comparing the outcome of the standard preprocessing pipeline with the outcome after incorporating prewhitening. Results from previous studies were replicated with the standard preprocessing pipeline but not with ones that included prewhitening. Importantly, the expected anticorrelation and antiphase state between ${\mathrm{HbO}}_2$ and HbR disappeared after prewhitening. Also, different prewhitening procedures yielded different results, both compared to the standard RS-NIRS preprocessing method and across prewhitening procedures.

These results open a discussion about which method would be more appropriate for those NIRS users who are interested in adding prewhitening to their preprocessing pipeline. Undoubtedly, the statistical challenges that have been recently described for both functional task based and connectivity NIRS data should be taken into consideration when analyzing NIRS data.^26,27 However, based on our results, we conclude that a better understanding of the effect of prewhitening in RS-NIRS data, and of the neurophysiological significance of the prewhitened signal, is still required to determine if prewhitening should be applied and, if so, which prewhitening procedure is more appropriate. Other relevant preprocessing methods might also be incorporated in the analysis of RS-NIRS data, such as global signal regression to attenuate the contribution of physiological confounds^18,23,56 or robust regression analysis to account for the presence of outliers due to motion artifacts in the data.^26,27 A full evaluation of the combined effect of these methods was beyond the scope of this paper, but we acknowledge the importance of including all these different methods in a potential discussion on standardized data preprocessing procedures for future RS-NIRS studies.