1.

Introduction

Functional near-infrared spectroscopy (fNIRS) is an emerging optical brain imaging technique that is becoming an increasingly popular method for pediatric brain imaging.¹ fNIRS measures the hemodynamic changes that effectively reflect brain activities occurring while people perform a wide range of mental tasks;^2–6 it can provide both topographic^3,4,7 and tomographic brain images.^6,8 Specifically, fNIRS monitors concentration variations of oxyhemoglobin (HbO) and deoxyhemoglobin (HbR) in the blood of the brain, through the skull, with absorption changes of near-infrared light at wavelengths between 700 and 1000 nm.⁹

Similar to electroencephalography (EEG), fNIRS allows for noninvasive observations of awake infants and children, but unlike EEG, fNIRS provides complimentary hemodynamic response data. Compared with functional magnetic resonance imaging, fNIRS devices are usually portable, silent, and use near-infrared light to detect brain activity.² These aspects of fNIRS technology make it a useful tool for researchers interested in child^10–20 and infant^21–27 development. Therefore, the use of this method with pediatric populations has grown threefold in the past decade.

Motion artifacts caused by participants’ head, jaw, eyebrow, or body movement can distort the data analysis.^28,29 A major impediment of interpreting fNIRS brain imaging data of child participants is that children’s data typically have more motion artifacts than the fNIRS data of adults. The field has put forth a variety of methods for addressing the fNIRS motion artifact problem in general.^28–30 The most direct method is trial rejection (for event-related design) or block rejection (for block design), which excludes all the data sections that contain motion artifacts from further analysis.³¹ Selb et al. compared several motion correction techniques on synthetic data and found trial rejection was the best method for preserving the original shape of the data. However, in the context of pediatric brain imaging, the major problem with this method is that the entire data samples tend to be relatively small, as children and infants have a short attention span and experimental designs can be as short as a few minutes. Thus, researchers have attempted methods that allow them to retain more of the collected data. These methods can be divided into two categories: the first category includes hardware-based methods (which need additional sensor or special optode setup), including short separation channel-based independent component analysis,³² using specially designed collodion-fixed fibers,³¹ accelerometer-based regression, or rescaling.^33,34 The difficulty with these methods is that additional hardware may further complicate and increase participant setup time, which is always a problem with young participants. The second category includes software-based methods, such as principal component analysis (PCA),³⁵ spline interpolation,³⁶ wavelet,³⁷ correlation-based signal improvement (CBSI),³⁸ and moving average (MA), plus Kalman filter and recursive least-square based methods for real-time preprocessing.^34,39

There have been a few prior studies that have attempted to explore the relative efficacy of different artifact correction approaches using simulated adult fNIRS data. These studies concluded that the wavelet method was best suited for the simulated adult fNIRS data.^28,29 Yet, it remains largely unknown whether the same is true of real pediatric fNIRS data from children subjects. Thus, in the present study, we aimed to improve the field’s methodology for child fNIRS brain imaging by directly comparing the efficacy of six software-based offline motion correction methods on real fNIRS data acquired from children (ages 6 to 12) completing a language task.

2.

Methods

2.3.

Functional Near-Infrared Spectroscopy Recordings

We used a TechEN-CW6 system with 690 and 830 nm wavelengths. The setup included one emitter and three detectors spaced 2.7 cm apart, yielding three data channels sampled at 10 Hz (Fig. 1). We examined brain activation in the left inferior frontal gyri (a language-related area). The probe localization was established and applied consistently for each participant using the international 10–10 transcranial system positioning;⁴³ Fz, Cz, and preauricular were measured for each participant and the emitter was anchored at F7 (see Fig. 1). The probes were mounted on the participants’ head using costume made caps (foam). Two caps were built to fit the different head sizes of the age range. Each cap was equipped with optode holders into which the optodes were plugged in. An additional wrapping band was used to secure the optodes in place.

Fig. 1

Measurement location and channel distribution.

2.4.

Data Processing

Data processing was completed using Homer2 fNIRS processing package⁴⁴ and an fNIRS data analysis tool (fNIRSDAT, a homemade software for general linear model (GLM) regression-based individual and group-level analysis) based on MATLAB^® (MathWorks, Natick, Massachusetts). In order to understand the performance of different methods, we categorized the motion artifacts in our data into four types: type A—defined as a spike with a standard deviation of 50 that was over the mean within a second; type B—a peak with a standard deviation of 100 from the mean during a time portion ranging from 1 to 5 s; type C—a gentle slope between 5 and 30 s with a standard deviation of 300 from the mean; type D—a slow baseline shifting longer than 30 s with a standard deviation of 500 from the mean.

The raw data were first converted into optical density change (delta OD) using the optical density change calculation function (hmrIntensity2OD in Homer2 software package) in Homer2. We then slightly modified the motion artifact identification function (hmrMotionArtifactbyChannel in Homer2 software package) and applied it to the delta OD data to detect the number of the type A motion artifact. The parameters were selected as $\mathrm{SDThresh}=50$, $\text{tMotion}=1$ (SDThresh, standard deviation threshold; tMotion, motion artifact time period); these parameters guaranteed the identification of the type A motion artifact (50 standard deviations to the mean within 1 s). However, the motion artifact identification function in Homer2 is only designed for identification of type A motion artifact (sudden move generated motion artifact). In this study, we also wanted to study the effect of type B, C, and D correction. Therefore, we used a home-made trend detection and detrending function using an MA algorithm to track the trend in fNIRS signal at different frequencies to identify type B, C, and D motion artifacts. This process also located the motion artifacts in the time series, thus providing the prior information for PCA and spline techniques.

We then applied five different motion correction techniques (PCA, spline, MA, wavelet, as well as MA and wavelet) to the delta OD series after the first round of the motion artifact detection process (except CBSI method; this method needs HbO and HbR data for motion correction; therefore, a different processing stream was applied for the CBSI method). Then, the motion artifact detection process was applied again in order to count the number of motion artifacts of different types. The motion corrected delta OD data were further filtered with a third-order Butterworth low-pass filter of $\sim 0.5\mathrm{Hz}$ to prevent the physiological noise interference (e.g., cardiac noises) and high-frequency measurement noises; the filter coefficients are automatically designed by the MakeFilter function in MATLAB^® (coefficients: $A$: [1 $-2.37$ 1.93 $-0.53$] B: [0.0029 0.0087 0.0087 0.0029]). Then the filtered OD data were converted to HbO and HbR concentration data by the modified Beer-Lambert law (MBLL) (differential path length factor is selected to be 6 for both 690 and 830 nm wavelengths). Finally, the HbO and HbR data derived from previous steps went through a block averaging process and were regressed with a GLM.⁴⁵ A flowchart of the processing steps can be found in Fig. 2.

Fig. 2

Data processing flowchart.

For CBSI method, we first applied a low-pass filter of ($\sim 0.5\mathrm{Hz}$) followed by the MBLL conversion to the delta OD data after the first round of motion detection process. The derived HbO and HbR data were then passed through the CBSI motion correction, and the processed HbO and HbR data were sent to the block averaging process and GLM regression, as well as converted back to delta OD data for the motion artifacts counting process (see Fig. 2).

3.

Results

The quantity variations of different motion artifacts are displayed in Fig. 3.

Fig. 3

The amount of variation of four types of motion artifacts after applying motion correction techniques.

For type A motion artifacts, the motion correction methods we applied, which produced the best results for both 690 and 830 nm wavelengths, were wavelet, CBSI, and the wavelet and MA combination, while the MA and PCA methods produced poor results for both wavelengths. For type B motion artifacts, the wavelet method alone, as well as the wavelet and MA combination, attained the best results for both wavelengths. For type C motion artifacts, the spline and wavelet and MA combination each performed better than other methods for both wavelengths. For type D motion artifacts, spline, MA, and the wavelet and MA combination outperformed the other methods for the 690 nm data, while all methods achieved the same performance compared to each other in 830 nm data, yet wavelet and CBSI performed better.

The GLM-based regression results are presented in Fig. 4. The $t$ values (averaged across three channels) of the HbO signal from all channels increased after we applied CBSI, wavelet, MA, and the wavelet and MA combination, but decreased after we applied PCA and spline. The $t$ values of the HbR signal changed in accordance with the HbO signal, but with a weaker magnitude. The $R$ values of both HbO and HbR decreased after we applied wavelet, MA, CBSI, and the wavelet and MA combination.

Fig. 4

General linear model regression derived $t$ value and $R$ value of three channels (three channels regressed as one) before and after applying different motion correction techniques.

Figure 5 shows the ${\mathrm{AUC}}_{0-1.5}$ results of the HbO and HbR data. All of the methods we applied achieved similar performances, except the spline and wavelet methods. Wavelet performed better than the other methods with a 65% decrease in HbO data and 64% increase in HbR data. Spline presented a less promising performance, showing only a 46% decrease in HbO data and a 46% increase in HbR data. The ${\mathrm{AUC}}_{\text{ratio}}$ results are shown in Fig. 6. The MA method achieved the best performance for both HbO and HbR data, while the spline and CBSI methods performed worse than other methods for both HbO and HbR data.

Fig. 5

${\mathrm{AUC}}_{0-1.5}$ variation after applying different motion correction techniques.

Fig. 6

${\mathrm{AUC}}_{\text{ratio}}$ before and after applying different motion correction techniques.

4.

Discussion

The goal of the study was to examine the efficiency of standard motion correction methods toward developmental fNIRS data. Specifically, we compared the performances of six motion correction techniques on fNIRS data derived from children 6 to 12 years old. The motion correction techniques we selected for comparison did not require extra equipment or a real-time data acquisition environment. Importantly, we did not consider the trial rejection method. The comparison metrics we selected aimed to qualitatively and quantitatively compare the performances of the motion correction abilities of these different techniques.

Previous studies reported that spline interpolation can achieve promising^28,36,47 and less promising results.²⁹ In theory, the spline technique can reduce all four types of motion artifacts equally. However, in our study, using this technique yielded minimal improvement in the data, in comparison to the data before motion correction techniques were applied in all metrics. The reason for this is that the spline technique uses the preidentified motion artifacts to subtract a generated cubic function. However, the motion artifact shapes vary in fNIRS data, thus causing unstable results. Vinette et al.⁴⁷ tested the spline interpolation method on real long-term fNIRS data measured from patients with epilepsy. Their results showed spline interpolation method is good at reducing baseline shift motion artifacts. However, they also pointed out the limitations of the spline interpolation method, requiring different parameters for different participants and lower performance when brain response signal and motion artifact overlapped. Thus, in our study, the use of same parameters for all participants and possible overlapping between task evoked responses and motion artifacts may lead to the less promising motion correction results.

Cooper et al. showed that using PCA produced better results than using no motion correction technique,²⁸ and Brigadoi et al. claimed the performance of PCA is variable depending on the data set.²⁹ In the current study, we examined a data set acquired from children. The performance of the PCA technique in our study was not promising. Our data show that PCA successfully reduced type B, C, and D motion artifacts, but did not reduce type A motion artifacts. Although the PCA technique performed fairly well in the area under the curve (AUC)-related metrics, the GLM-based regression-related metrics decreased after applying PCA for both HbO and HbR. The reason for this may be that the PCA technique is a multichannel approach; thus, it requires that an artifact be presented in multiple channels (to be identified as a principal component), which is not always the case. The basic assumption of the PCA technique, that the components need to be independent, was violated. This may have additionally contributed to the poor performance we observed.

The CBSI technique performed well with our data. It decreased the number of all motion artifacts, except type B. The ${\mathrm{AUC}}_{0-1.5}$ and ${\mathrm{AUC}}_{\text{ratio}}$ metrics jointly showed an improvement in both averaged HbO and HbR data. The GLM-related metrics increased in $t$ value, and also slightly decreased in residual for both HbO and HbR data. However, the results of CBSI can be unstable, especially because an assumption of the CBSI method is that the HbO and HbR data are always negatively correlated, which may not always be the case. According to our data, in most cases, when the HbO signal increases, the HbR signal decreases, but sometimes HbR increases or remains unchanged. Furthermore, this is a reason why fNIRS researchers have reported that the HbO signal is more reliable than the HbR signal in fNIRS studies.⁴⁸

In our study, the wavelet technique achieved a promising result, as it greatly reduced type A and type B motion artifacts, and barely reduced type C and type D motion artifacts. It also achieved very good results in two AUC-related metrics, and slightly increased the $t$ value and $R$ value for both HbO and HbR data. An advantage of the wavelet technique is that it does not require any prior motion artifact detection algorithm. Wavelet decomposes the signal into components on different time resolutions and then detects and fixes motion artifacts, thus eliminating the need for detection algorithms. Finally, the corrected component signals are restored to form the motion corrected signal. Brigadoi et al.²⁹ concluded that the wavelet technique achieved the best results in their study. However, according to our data, wavelet showed a weaker ability in reducing type C and type D motion artifacts. Therefore, we introduced the MA technique to work jointly with the wavelet technique to enhance the motion correction results.

Our data showed that using MA accomplished decent performance on reducing type C and type D motion artifacts. It also achieved good results for all other four metrics. The problem with the MA method is that it has a weak ability to reduce type A and type B motion artifacts. Therefore, the MA method may not work independently when type A and type B motion artifacts are significantly dominant in the data.

The wavelet and MA methods independently yielded some of the best results in reducing type A/B and type C/D motion artifacts. Due to this, we applied wavelet and MA combined methods to our data to observe the cumulative results. As expected, the wavelet and MA combination method attained the best performance of reducing motion artifacts in our data. The method greatly reduced all four kinds of motion artifacts and acquired very good AUC-related metrics. The GLM regression $t$ value for HbR signal we obtained was not as high as CBSI method. The $R$ value we obtained was the smallest among all methods, indicating that this motion corrected data had the smallest deviation from the modeled data. As stated above, the reason for this is that the wavelet method and the MA method can, respectively, handle high-frequency and low-frequency noises. However, the disadvantage of using these combined methods is that the motion correction process must be applied two times, which may cause an over-reduction of the magnitude of the signal and can lead to low regression coefficients and low $t$ values.

In general, the spline technique was fairly successful at reducing all types of motion artifacts, but leads to unsatisfied metrics. This was possibly due to improper correction of motion artifacts. The PCA method reduced different motion artifacts less successfully than other methods in this study, thus leading to less promising results in all metrics except ${\mathrm{AUC}}_{0-1.5}$. The reason for this could be the absence of motion artifacts in the channels in our study. If so, PCA could not catch the motion artifacts as a principal component. CBSI demonstrated good results for reducing type A and type C motion artifacts, wavelet reduced type A and type B motion artifacts better than other methods, MA corrected type C and type D motion artifacts better than other methods, and the combination of wavelet and MA effectively reduced all types of motion artifacts. These latter four techniques achieved similar results in ${\mathrm{AUC}}_{0-1.5}$, indicating they reduced the curve from 0 to 1.5 s. MA and wavelet combination and MA performed better than other methods for the ${\mathrm{AUC}}_{\text{ratio}}$ metric. Interestingly, both MA and wavelet combination and MA achieved promising results for $t$ value and $R$ value metrics. This indicates that type C and type D motion artifacts affect GLM metrics more than type A and type B motion artifacts. Table 2 provides a performance summarization for all the methods applied in the current study.

Table 2

Methods’ performance summarization.

5.

Conclusion

In this study, we applied six different techniques on fNIRS data recorded from children completing a language task and compared the performance of these techniques. We compared the differences in amount of motion artifacts, ${\mathrm{AUC}}_{0-1.5}$, ${\mathrm{AUC}}_{\text{ratio}}$, GLM regression-based $t$ values and $R$ values before and after motion correction application. Our data revealed that the MA and wavelet combination produced the largest improvement in the computed metrics. We also found that type C and type D motion artifacts affect GLM-related metrics more than type A and type B motion artifacts.