2.1.

Subjects and Data Collection

NIRS data were collected in two distinct groups: 46 healthy adults and 36 stroke patients. Forty-six healthy volunteers with no history of cardiovascular or neurological disease were recruited at the University of Copenhagen, Denmark, and 36 patients with symptoms of stroke in the territory of the internal carotid artery were recruited from the Stroke Unit, Glostrup Hospital in Denmark. The Biomedical Research Ethics committee in the Capital Region of Denmark approved both studies (H-B-2008-039 and H-B-2008-088). All subjects gave informed consent to participate in the study, which was undertaken in accordance with the Helsinki Declaration of 1964, as revised in Edinburgh in 2000. Subsets of this data have been previously analyzed and published independently.^21,30,38 The NIRS data were collected with a TechEn Inc. (Milford, Massachusetts) CW6 system. NIRS optodes were placed bilaterally and symmetrically on the forehead, with, on each side, one source (690 and 830 nm) and one detector 3 cm away, avoiding the midline sinus. Additional channels were recorded at other locations in the two groups, but were not analyzed for the present study. NIRS recordings were acquired at a sampling rate of 200 Hz and were subsequently downsampled to 25 Hz.

Each dataset consisted of 10 min of resting-state NIRS recording. Both the healthy subjects and stroke patients were placed in a supine position in a silent room with a constant temperature and the light dimmed. After 15 min of rest in the supine position, a 10-min trial of spontaneous breathing was acquired.

2.2.

Data Analysis: Interhemispheric Correlation

For each channel, the raw intensities at 690 and 830 nm were converted to relative changes in optical density $\mathrm{\Delta }\mathrm{OD}$ using the modified Beer–Lambert law with a differential path length factor of 6 at both wavelengths. When applicable, MAs were identified and corrected on the $\mathrm{\Delta }\mathrm{OD}$ time traces as detailed in Sec. 2.3. Then the corrected $\mathrm{\Delta }\mathrm{OD}$ time series were converted to changes in oxy and deoxyhemoglobin concentrations, $\mathrm{\Delta }[\mathrm{HbO}]$ and $\mathrm{\Delta }[\mathrm{HbR}]$, respectively. The hemoglobin concentration time series were band-pass filtered (third-order Butterworth low-pass filter and fifth-order Butterworth high-pass filter) in four different frequency bands of physiological interest: cardiac (0.5 to 1.5 Hz), respiration (0.15 to 0.5 Hz), low-frequency oscillations (LFOs, 0.05 to 0.15 Hz), and very LFOs (VLFO, 0.01 to 0.07 Hz). Finally, our metric of interest was the IHC between the two symmetrical channels on the forehead. Specifically, the IHC was computed as the correlation coefficient over the whole 10 min recording duration between the filtered $\mathrm{\Delta }[\mathrm{HbO}]$ time series of the two symmetrical channels on the forehead for each of the four frequency bands. Therefore, for each subject’s 10 min recording, we obtained four IHC values corresponding to the four physiological frequencies. The IHC values are between $-1$ and $+1$, with $+1$ meaning that the symmetrical channels are perfectly in phase, $-1$ representing perfect anticorrelation (180-deg phase), and 0 representing either a complete lack of correlation between the channels or a 90-deg phase shift.

2.3.

Motion Artifact Identification and Correction

MAs were detected on the $\mathrm{\Delta }\mathrm{OD}$ time series using the function hmrMotionArtifactPerChannel from the Homer2 processing package,^39,40 except for the wavelet motion correction procedure, which performs the MA identification as part of the correction analysis. For each data channel, the algorithm automatically identifies MAs based on fast changes in signal amplitude and/or standard deviation. If, over a temporal window of length tMotion, the standard deviation increases by a factor exceeding SDThresh, or the peak-to-peak amplitude exceeds AMPThresh, then the segment of data of length tMask starting at the beginning of that window is defined as motion. We used $\mathit{\text{tMotion}}=0.5\mathrm{s}$, $\mathit{\text{tMask}}=2\mathrm{s}$, $\mathit{\text{SDThresh}}=14$, and $\mathit{\text{AMPThresh}}=0.2$ for all datasets. In practice, these parameters may be tweaked on an individual-subject or -study basis in order to optimize the MA identification as confirmed by visual inspection. In the present study, in order to perform an objective characterization of MA incidence and of their correction, we fixed these parameters at the same values for all datasets in both groups.

We compared five different MA correction approaches, which are illustrated in Fig. 1, left: not correcting the MAs (“No correction”), discarding the MA segments (“Discard”), spline interpolation correction (“Spline”),¹⁶ recursive-targeted PCA correction (“tPCA”),¹⁵ and wavelet filtering (“Wavelet”).¹⁸ The last three MA correction methods have been fully described by their respective authors and have been summarized in previous publications,^15,21,22 so we only give a brief overview of the different algorithms here.

Fig. 1

Example illustrating: (a) a motion-free segment of data, (b) the same data after addition of MA segments, which are identified by the gray areas, and the different motion correction procedures: (c) discarding the MA segments, (d) spline interpolation, (e) targeted PCA, and (f) wavelet filtering. Left: $\mathrm{\Delta }\mathrm{OD}$ at 830 nm for the symmetrical forehead channels. Right: $\mathrm{\Delta }[\mathrm{HbO}]$ time series for the symmetrical channels, after filtering in the LFO frequency band (0.05 to 0.15 Hz). The displayed IHC was computed over the whole 600 s of recording, but only 160 s of data are shown for better visualization of the MAs.

2.4.

Motion Artifact Addition

A subset of the healthy datasets was used as reference data with minimal motion contamination. We selected all healthy datasets for which the cumulative duration of MA segments was less than 5% of the total recording time (34 datasets). On these “motion-free” datasets, we added MA segments obtained from the stroke datasets.

2.4.1.

Same motion artifact content as the stroke group

For each healthy recording, we randomly chose 20 stroke recordings from which to extract the MA segments, yielding a total of 680 datasets with added MAs. For each recording, the procedure to add MA segments, illustrated in detail in Fig. 2, was the following: one clean healthy dataset with no MA [Fig. 2(a)] was matched with a MA-contaminated stroke dataset [Fig. 2(b)]. The MA segments were identified on the stroke time series [Fig. 2(b)] as described in Sec. 2.3. Then the MA segments from the stroke $\mathrm{\Delta }\mathrm{OD}$ time series were substituted in the healthy $\mathrm{\Delta }\mathrm{OD}$ [Fig. 2(c)] on a per-channel and per-wavelength basis. Importantly, each MA segment was shifted in the healthy data by the same offset it was shifted in the stroke data. This shift can, for instance, be observed around $t=120\mathrm{s}$ on Figs. 2(b) and 2(c). An additional example of the procedure was provided in Figs. 1(a) and 1(b). The procedure was then repeated for all 20 selected stroke recordings.

Fig. 2

Illustration of the addition of motion artifacts (MAs) to healthy datasets: (a) The original $\mathrm{\Delta }\mathrm{OD}$ time series from a healthy subject, with no identified MA. (b) The MA segments are identified on the $\mathrm{\Delta }\mathrm{OD}$ time series of a stroke patient. (c) The MA segments are substituted in the healthy time series. The process is performed per source-detector pair and per wavelength. Only one channel is shown here (left side, 830 nm). The data offsets resulting from the MAs are also incorporated in the modified data. See for instance around $t=120\mathrm{s}$, the baseline $\mathrm{\Delta }\mathrm{OD}$ has been shifted negatively by about 0.2 after the MA segment to reflect the actual effect of the motion observed in the stroke dataset.

2.4.2.

Larger motion artifact content than the stroke group

In addition to the previous procedure where we artificially introduced the same MA content as that of the stroke datasets, we studied the effect of introducing further MA contamination. For this procedure, we first created a “pool” of 355 MA segments extracted from the stroke datasets as described above. Each MA consisted of the four segments of $\mathrm{\Delta }\mathrm{OD}$ at the two wavelengths and the two channels. For each healthy dataset, we then added MA segments randomly selected from this pool at randomly selected time points, following the same segment substitution as described above (Sec. 2.4.1), in order to obtain a target average MA contamination. For each healthy dataset, we created 20 instances of randomly placed MA segments for the target MA contamination level. Finally, the whole process was repeated for different MA contamination levels varying from 17% up to 50% (duration of the cumulated MA segments over total recording duration). Figure 3 shows an example of a $\mathrm{\Delta }\mathrm{OD}$ time series from the healthy dataset after introduction of a motion segment that results in approximately 10%, 30%, and 50% MA content. Note that this procedure results in a large number of short MAs as opposed to few long segments of MAs, so that the contamination higher than 50% corresponds to unusable datasets.

Fig. 3

Illustration of the addition of different levels of motion contamination to healthy data. The original $\mathrm{\Delta }\mathrm{OD}$ time series from a healthy subject (blue) is contaminated by substitution of MA segments obtained from the stroke datasets, resulting in cumulative duration of MAs of 10%, 29%, and 52% (cyan).

2.5.

Statistical Analysis

The distribution of the number of MA segments is not Gaussian, and the distribution of IHC is bounded between $-1$ and $+1$ and is non-Gaussian. We, therefore, report the median and IQR of the number of MA segments and of the IHC values in the four frequency bands of interest for each group.

To quantify the dependence of IHC on the number of MAs, we computed the Pearson correlation coefficient between the number of MAs and the Fisher’s Z-transform of the IHC. The Fisher Z-transformation normalizes the distribution of the IHC.

We compared the IHC values between the following datasets: healthy without MA versus healthy with added MAs; healthy without MA versus healthy with added MAs, after MA correction; stroke versus healthy without MA; and stroke versus healthy with added MAs, before and after MA correction in both groups. To compare the IHC values between two datasets, we performed a two-sample Student $t$-test applied on the Fisher Z-transformation of the IHC. Note that we report the median and IQR of the IHC before transformation for convenience, but all statistical analyses were performed on the Fisher Z-transform of the IHC. The threshold for significance was set to $p<0.01$.

2.6.

Overall Method

In summary, our approach was the following. First, in Sec. 3.1, we study the correlation between incidence of MAs, and IHC, on the original (i.e., not MA-corrected) healthy and stroke datasets. Second, in Sec. 3.2, using only the healthy datasets with minimal MA contamination, we investigate the effect on IHC of adding MAs and of correcting them with the different approaches. We study in detail the case where synthetic motion contamination is equivalent to that of the stroke group, and further investigate how the different MA correction approaches perform when the MA contamination increases to up to 50% of the total recording time. This section provides a case where the ground truth IHC is known, and we can directly study the effect of MAs and their correction. Finally, in Sec. 3.3, we compare the healthy group after the addition of MA and the stroke group, in both cases with or without MA correction. The objective of this last part is to compare the healthy and stroke groups with similar MA contamination so as to minimize the confounding effect of motion contamination.

3.1.

Interhemispheric Correlation Correlates with the Number of MAs

As illustrated in Fig. 4(a), without MA correction, the IHC (median and IQR) was significantly higher in the healthy group than in the stroke group at all frequencies: cardiac (healthy: 0.98, 0.95 to 0.99; stroke: 0.93, 0.84 to 0.96; $p=5\times {10}^{-6}$), respiration (healthy: 0.78, 0.55 to 0.90; stroke: 0.55, 0.22 to 0.76; $p=3\times {10}^{-5}$), low-frequency (healthy: 0.78, 0.66 to 0.89; stroke: 0.52, 0.01 to 0.68; $p=2\times {10}^{-7}$), and very-low frequency (healthy: 0.78, 0.64 to 0.90; stroke: 0.57, 0.41 to 0.71; $p=3\times {10}^{-5}$). The IQR values show that the variability across subjects was higher in the stroke group than in the healthy group in all frequency bands [see also error bars on Fig. 4(a)].

Fig. 4

(a) IHC in the healthy group and the stroke group obtained with the original datasets (i.e., without MA addition or correction), in the four frequency band of interest. The colored circle represents the median and the error bars the IQR over all subjects. (b) Number of MA segments in the healthy group and the stroke group, obtained from the original datasets (median and IQR). For (a) and (b), the significance of the difference between the two group is displayed as *$p<0.01$ and **$p<0.001$. (c) Scatter plot of the Fisher’s Z-transform of the IHC versus the number of MAs for the healthy subjects and the stroke patients. The displayed Pearson correlation coefficient $R$ and corresponding $p$-value were computed over the combined data (healthy and stroke).

The number of MA segments was significantly lower ($p=0.0019$) in the healthy group ($\text{median}=2$ and $\mathrm{IQR}=1$ to 5) than in the stroke group ($\text{median}=8$ and $\mathrm{IQR}=4$ to 13), as displayed in Fig. 4(b).

Figure 4(c) shows the scatter plot, combined for both healthy and stroke groups, of the IHC in the low-frequency band versus the number of MAs. There is a significant negative correlation between the number of MAs and the IHC value (Pearson correlation coefficient: $R=-0.38$, $p=5\times {10}^{-4}$). The correlation was also significant in the cardiac ($R=-0.42$, $p=1\times {10}^{-4}$), and respiration ($R=-0.30$, $p=6\times {10}^{-3}$) frequency bands, and marginally significant in the VLFO ($R=-0.27$, $p=0.015$) frequency band.

3.2.

Adding Motion Artifacts to Motion Artifact-Free Datasets Significantly Reduces Interhemispheric Correlation

Thirty-four healthy datasets were considered motion-free as defined by the total duration of MA segments being under 5% of the total recording time (i.e., under 30 s over 10 min of recording). On this subset of data, the median (IQR) IHC was 0.98 (0.95 to 0.99), 0.85 (0.65 to 0.91), 0.80 (0.71 to 0.91), and 0.81 (0.66 to 0.90) in the cardiac, respiratory, LFO, and VLFO frequency bands, respectively. As expected from the results above, these values are slightly higher than the median values for the complete healthy group.

After the addition of MAs on the healthy data subset, the IHC decreased significantly to 0.88 (0.77 to 0.95), 0.49 (0.23 to 0.73), 0.56 (0.34 to 0.74), and 0.60 (0.35 to 0.79) at the cardiac, respiration, LFO, and VLFO frequencies, respectively. This is significantly lower than the original values, with a $p$-value of $1\times {10}^{-6}$, $3\times {10}^{-12}$, $7\times {10}^{-9}$, and $8\times {10}^{-7}$, respectively. This result is illustrated in Fig. 5 by the dark blue and cyan bars. As can be observed from the IQR values and from the error bars on Fig. 5, the addition of MAs also increased the variability of the IHC across subjects.

Fig. 5

Bar graph of the IHC in the healthy datasets with less than 5% MA contamination (blue), after addition of MA segments (cyan), and after motion correction using the discard method (gray), spline interpolation (green), targeted PCA (yellow), and wavelet analysis (magenta). The bar height represents the median, and the error bar the IQR over all subjects.

3.3.

Effects of Different Motion Correction Approaches

We performed a sensitivity analysis on the input parameters $p_{\text{Spline}}$, ${\sigma }_{\mathrm{PCA}}$, and ${\alpha }_{\mathrm{IQR}}$ for the Spline, tPCA, and Wavelet approaches, respectively. For each correction method, we assigned three values to the input parameters as described in Sec. 2. Working with the healthy dataset, we selected the value of each parameter that resulted in the least difference between the original IHC (on the MA-free datasets) and the IHC after MA addition and correction. This resulted in $p_{\text{Spline}}=0.99$, ${\sigma }_{\mathrm{PCA}}=97\%$, and ${\alpha }_{\mathrm{IQR}}=1.5$. These values were employed in all subsequent results below.

Figure 1 shows an example of the effect of the MA addition and correction on a short segment of motion-free healthy data. The first column displays the $\mathrm{\Delta }\mathrm{OD}$ time traces at 830 nm of the two symmetrical channels, while the second column displays the corresponding $\mathrm{\Delta }[\mathrm{HbO}]$ times series after filtering in the 0.05 to 0.15 Hz frequency band.

Figure 5 shows the effect of the different motion correction methods in the four frequency bands of interest for the healthy datasets with added MAs. All approaches increased the median IHC compared to the dataset with added MAs. However, all of them retrieved a value lower than the true IHC. Indeed, the IHC remained significantly lower than their original true values for all correction approaches.

The performance of the different approaches differed. The discard approach performed best in all frequency bands, increasing the IHC to 0.95 (0.91 to 0.98) in the cardiac frequency band, 0.72 (0.50 to 0.85) in the respiration frequency band, 0.71 (0.56 to 0.81) in the LFO band, and 0.69 (0.51 to 0.79) in the VLFO band. The spline method increased the IHC to 0.91 (0.84 to 0.95) in the cardiac frequency band, 0.62 (0.38 to 0.78) in the respiration frequency band, 0.66 (0.48 to 0.77) in the LFO band, and 0.66 (0.50 to 0.78) in the VLFO band. Targeted PCA retrieved an IHC of 0.91 (0.82 to 0.96) in the cardiac frequency band, 0.57 (0.32 to 0.78) in the respiration frequency band, 0.61 (0.38 to 0.78) in the LFO band, and 0.63 (0.41 to 0.80) in the VLFO band. Finally, the wavelet analysis increased the IHC to 0.92 (0.86 to 0.96) in the cardiac frequency band, 0.71 (0.52 to 0.82) in the respiration frequency band, 0.65 (0.49 to 0.75) in the LFO band, and 0.62 (0.40 to 0.79) in the VLFO band.

3.4.

Effect of Higher Motion Artifact Content

Figure 6 shows the evolution of the IHC, without MA correction and for the different MA correction approaches, as a function of the added MA content. The IHC decreases with increased MA content up to approximately 45% MA contamination. For higher MA content, the IHC starts to increase sharply. All MA correction methods improve the IHC, with the exception of the tPCA approach which degrades the IHC for MA content higher than $\sim 25\%$. The discard approach performs best at all MA contents, and wavelet and spline have similar performances. However, even after MA correction, the IHC decreases for increasing MA contamination.

Fig. 6

IHC in the healthy datasets with less than 5% MA contamination after added MA, versus level of added MA contamination (expressed as percent of the total recording time). The results are presented after addition of MA segments (cyan), and after motion correction using the discard method (gray), spline interpolation (green), targeted PCA (yellow), and wavelet analysis (magenta). The blue and red arrows indicate the level of MA contamination corresponding to the original healthy datasets (no MA added) and to the original stroke datasets.

3.5.

Comparison of Healthy with Added Motion Artifacts, and Stroke

Figure 7 shows the effect of the different motion correction methods in the four frequency bands of interest for the stroke datasets. After the Discard motion correction, the median (IQR) IHC in the stroke group increased to 0.97 (0.93 to 0.98) at the cardiac frequency, 0.70 (0.51 to 0.77) at the respiration frequency, and decreased to 0.48 (0.34 to 0.66) in the LFO frequency band, and to 0.50 (0.34 to 0.62) in the VLFO frequency band. The spline procedure increased the IHC to 0.94 (0.88 to 0.96) at the cardiac frequency and 0.57 (0.42 to 0.71) at the respiration frequency, and decreased it to 0.46 (0.27 to 0.64) in the LFO frequency band, and 0.47 (0.29 to 0.61) in the VLFO frequency band. tPCA increased the IHC to 0.94 (0.82 to 0.97) at the cardiac frequency, and 0.62 (0.24 to 0.75) at the respiration frequency, and decreased it to 0.49 (0.27 to 0.67) in the LFO frequency band, and 0.55 (0.47 to 0.75) in the VLFO frequency band. Finally, the wavelet filtering increased the IHC to 0.96 (0.91 to 0.97) at the cardiac frequency, and 0.67 (0.47 to 0.76) at the respiration frequency, and had little effect in the lower frequency bands, where it retrieved 0.50 (0.30 to 0.67) in the LFO frequency band, and 0.57 (0.42 to 0.72) in the VLFO frequency band.

Fig. 7

Bar graph of the IHC in the stroke datasets (red), and after motion correction using the discard method (gray), spline interpolation (green), targeted PCA (yellow), and wavelet analysis (magenta). The bar height represents the median and the error bars the IQR over all subjects.

The IHC after MA addition in the healthy dataset was not significantly different than the IHC in the stroke group, in the respiration ($p=0.42$), LFO ($p=0.08$), and VLFO ($p=0.59$) frequency bands. At the cardiac frequency, the IHC was marginally lower in the healthy dataset with added MA than in the stroke dataset ($p=0.02$). These results are illustrated in Fig. 8.

Fig. 8

Comparison of the IHC in the healthy group and the stroke group. The original healthy datasets (i.e., before MA addition, in dark blue) present high IHC in all frequency bands. After addition of MA segments obtained from the stroke dataset, the IHC in the healthy group (cyan) is significantly reduced at all frequencies and is not significantly different from the stroke IHC (red). After applying the discard procedure to these datasets (healthy with added MAs and stroke), the IHC is significantly lower in the stroke group than in the healthy group in the LFO frequency band, and marginally lower in the VLFO frequency band. Significance of the difference between the healthy and stroke datasets: (o) not significant (i.e., $p>0.01$), *$p<0.01$, **$p<0.001$.

Figure 8 also illustrates the comparison between the healthy group (without MA and with added MAs), and the stroke group, before and after the discard correction. After the discard motion correction, the IHC in the stroke group and in the healthy group with added MAs were not significantly different in the cardiac ($p=0.37$) and respiration ($p=0.77$) frequency bands. The IHC was significantly lower in the stroke group for the LFO ($p=3\times {10}^{-5}$), and for the VLFO ($p=1\times {10}^{-4}$) frequency bands. We are not showing the comparison between the stroke and healthy groups for the other motion correction approaches, as they performed worse than discard in retrieving the “true” IHC in the healthy datasets, as detailed in Sec. 3.3.

4.1.

Negative Correlation of Interhemispheric Correlation and Number of Motion Artifacts

We observed that the IHC at all frequencies in the original datasets had a negative correlation with the number of MAs (results only shown for the LFO frequency band). Furthermore, adding MA segments to the healthy datasets with minimal MA contamination resulted in a significantly reduced IHC at all frequencies and in an increased IHC variability across subjects.

The reduction of IHC with increased MA content can be counter-intuitive. If we assume motion to induce similar artifacts on all channels, we would expect the correlation between symmetrical channels to increase artificially with the presence of MAs. This effect was sometimes observed (data not shown). However, in most cases, the fast signal changes due to MA occur in opposite directions in different channels or with a slight delay, inducing a decrease in the IHC, as can be clearly observed in Fig. 1(b). See, for instance, in Fig. 1(b), right, around $t=350\mathrm{s}$, an MA is present on both channels, resulting in an increase of the [HbO] fluctuations, albeit with a short delay between the two channels, which reduces IHC. The combination of these two effects results in an overall decrease of the IHC, but with increased variance. Note, however, that for very high MA content, the IHC starts to increase again, as illustrated in Fig. 7 for MA content $>45\%$. The very high number of synthetically incorporated MA segments is inducing correlations between channels. Due to the nature of the MA contamination we model here, which is typically a large number of short (a few seconds) motion segments, 50% MA contamination represents extremely “noisy” datasets (one out of two data points identified as motion-corrupted), which would be considered unusable in practice.

Similar observations have been made in rs-fc blood oxygenation level-dependent fMRI studies. Power et al.^41,42 reported that the motion affected correlation patterns over the brain, even when compensatory spatial registration and regression of motion estimates from the data were employed. They also observed that, while MAs tend to induce spurious increased correlations, some long-distance correlations were instead decreased by subject motion.⁴¹ Here, we observe an overall decrease in symmetrical forehead correlation with increased subject motion, but accompanied by increased variance, revealing that the IHC increased in some subjects and decreased in others.

The contamination by MAs is an important confounding factor to take into account when interpreting NIRS data analysis based on correlations, as is mostly the case for oscillation-type studies such as functional connectivity and autoregulation studies. It is particularly crucial to be aware of this effect when comparing different populations, or different paradigms, that are prone to different levels of motion contamination. In the study presented here, the analysis of the original datasets revealed a highly significantly lower IHC in the stroke group compared to the healthy group at all frequencies. However, the number of MAs was also significantly higher in the stroke group, due to the difficulty for this clinical population to keep their head still during the whole duration of the recording. This suggested that the lower IHC in the stroke group was at least partly due to the higher MA content. Indeed, when adding motion segments identified on the stroke group to the healthy datasets, the IHC decreased to the level of the stroke group.

In general, studies comparing a clinical population to a control healthy group might be subject to the same confounding effect of different motion contaminations. This possible confound should also be taken into account when studying the dependence of oscillations or functional connectivity on functional activity.⁴³ A functional paradigm involving a motor or speech task could induce more MAs than the resting condition. Based on the present work, we, therefore, recommend that every study using correlation-based metrics and making comparisons between subjects or over time should test for MA content.

4.2.

Implementation of the Discard Approach

The implementation of the discard approach requires some care in the data analysis streamline. This is because the band-pass filtering step can spread artifacts due to motion from MA-contaminated segments onto the adjacent MA-free segments. If the band-pass filtering was performed before the removal of the MA segments, the large signal changes due to MAs would spread into the adjacent good MA segments through the filtering process. However, if instead the MA segments were first removed and the remaining good segments concatenated, this would also artificially modify the frequency content of the time series. Similarly, if the MA-free segments were zero-padded, spurious frequency content would artificially arise at each transition between MA-free and MA-contaminated segments replaced with zeros. Instead we first identify the MA segments and correct them using the spline interpolation approach in order to minimize any large change in signal. We then apply the band-pass filtering, and finally compute the IHC only on the segments originally identified as motion free.

This same issue has already been recognized for the correction of MAs in rs-fc fMRI data and has been recently thoroughly discussed by Power et al.^41,42 They proposed the following approach for fMRI resting state data analysis: first the “bad” data points contaminated by MAs are identified using objective data quality metrics, then the data points of these MA segments (the censored “bad” time points) are reconstituted using the MA-free segments (the uncensored “good” data), using a least-square spectral decomposition adapted for nonuniformly sampled data.⁴⁴ Next, the reconstituted signal is band-pass filtered. Following the frequency-filtering, the reconstituted and filtered signals are “recensored,” i.e., the data points corresponding to the originally identified MA segments are ignored. Power et al.⁴² refer to this procedure as “scrubbing” the data. We implement a similar approach here, but we reconstitute the signal using the spline interpolation instead of the least square spectral decomposition. The rationale is identical: we implement a motion correction routine on the MA-contaminated segments prior to band-pass filtering, because, even though these segments are ultimately not included in the IHC computation, the frequency filtering step would otherwise spread the MA effect to adjacent good segments.

4.3.

Effect of Motion Artifact Correction in Different Frequency Bands

We observed varying performances of the different motion correction methods. In all frequency bands, the best approach was to discard the segments with MAs. This good performance when discarding the MA segments stands in contrast to results for evoked response to activation. Indeed, our previous publications^15,21,22 have shown that for evoked response studies, it is always preferable to correct the MAs than to reject the stimulus epochs contaminated by MAs. This may be simply explained by the amount of information ignored by discarding the MAs. For evoked response studies, the HRF is generally estimated, through block-averaging or deconvolution, over a relatively small number of stimulus epochs. When an MA is detected within a stimulus epoch, the whole block of data (typically 5 to 30 s or more) is rejected. In contrast, in the present study, only the data corresponding to the duration of the MA are rejected, typically several short segments of a few seconds which cumulative duration is on the order of 10% of the total recording.

Note that in the present results, we did not correct the retrieved IHC values for the number of time points included in the correlation computation, which is lower for the discard approach than for the other methods. We verified that this had little influence on the retrieved IHC for the level of MA contamination investigated here. Specifically, for all healthy datasets, we computed the IHC over (a) the whole duration of the recording (i.e., 10 min), and (b) including only subsets of time points corresponding to the MA-free segments in the stroke datasets. For each healthy dataset, we repeated the procedure for 20 randomly chosen stroke datasets. The retrieved IHC differed by less than 0.3% in all frequency bands when using all time points or only a subset. This gives us confidence that the number of time points has little influence on the retrieved IHC values for MA contamination on the order of 10%, typical of clinical data as reported in the present study.

Second to discarding the MA segments, the wavelet approach worked reasonably well for the high frequency bands (cardiac and respiration). However, the Wavelet approach performed poorly in the low and very low frequency bands. We believed that this is due to the fact that the wavelet approach does not correct sudden offsets in the NIRS signal. Instead the wavelet analysis acts as a low-pass filter so that the abrupt signal change is transformed into a slower drift. This drift has little impact on the signal at the cardiac and respiration frequencies, but introduces spurious frequency content in the LFO and VLFO bands.

Conversely, spline performed best in the low frequency bands (LFO and VLFO, second best approach after discard), but poorly at high frequencies (cardiac and respiration). From visual inspection of the data after correction [see for instance Fig. 1(d)], we believe that this is due to the fact that the spline interpolation models the signal “too well” over the MA segment, i.e., it also models the high frequency physiology, which is then subtracted from the segment. For each MAs segment, the cardiac and respiration physiologies are removed by the spline approach, thus reducing the IHC in these frequency bands. In contrast, this subtraction over a few seconds has little impact on the low frequency physiology.

Finally, tPCA performed relatively poorly in all frequency bands. This may be due to the limited number of channels (only two channels at two wavelengths each) in the present study. A future study should investigate whether tPCA performs better in oscillations studies with more channels.

The frequency bands of interest for functional connectivity and autoregulation studies are generally the LFO and VLFO bands. In these frequency bands, we suggest that correcting the MA should be achieved by discarding the MA segments following the approach described here, or as a second best approach with spline interpolation. We need to emphasize, however, that even after MA correction of the added MA, the IHC remains significantly smaller than the original values. This result suggests that IHC data should still be interpreted with caution after motion correction.

4.4.

Difference Between Healthy and Stroke

While the main objective of this study was to investigate the effect of MAs on resting-state data, we also observed an interesting physiological finding. The IHC was smaller in the stroke group than in the healthy group before MA correction. However, we showed that this lower IHC could at least be partly attributed to the higher MA content. After MA correction, this difference persisted but was less significant. It is unclear if the remaining difference was due to a real underlying physiological effect or was an artifact arising from suboptimal MA correction, or from possible other confounds.

Therefore, we compared instead the stroke group to the healthy group after MA addition, and corrected MA segments with the discard procedure in both cases. After this process, the IHC was not significantly different between the two groups at the cardiac frequency ($p=0.37$) and at the respiration frequency ($p=0.77$). The IHC was significantly lower in the stroke group in the LFO frequency band ($p=3\times {10}^{-5}$) and in the VLFO frequency band ($p=1\times {10}^{-4}$).

This finding suggests that the remaining difference arises from a real physiological effect, possibly originating from impaired autoregulation on one hemisphere, or from decreased functional connectivity across hemispheres. In agreement with this explanation is the fact that the difference in IHC is observed at the low and very low frequencies, around 0.1 and 0.04 Hz, but not at the higher cardiac and respiration frequencies, for which dynamic autoregulation is too slow to operate. Indeed, within the simplified framework of transfer function analysis, dynamic cerebral autoregulation is considered to act as a high pass filter, where fast changes in blood pressure are translated into fast changes in cerebral blood flow, while slower fluctuations are attenuated and delayed.^45,46 Further studies will be required to elucidate the exact mechanism, whether of neuronal or vascular origin, behind our finding.

4.5.

Limitations of the Study