2.1.

Sample and Data Acquisition

Eighteen healthy native French-speaking adults participated in this study. One participant was excluded from the analyses, as it showed continuous noise precluding the identification of individual artifacts. The final sample for the validation of using PARAFAC for artifact correction thus included 17 participants (mean age ± standard deviation = 22.8 ± 2.0 years; nine females and eight males). All were right-handed and presented no neurological or psychiatric disorders. Experimental procedures were approved by the local ethics committee.

fNIRS data were acquired with a multichannel Imagent Tissue Oxymeter (ISS Inc., Champaign, Illinois) frequency-domain fNIRS device using 14 light detectors and 60 laser light emitters, each regrouping two light sources of different wavelengths ($\lambda 1|2=690|830\mathrm{nm}$) with an average power of 10 mW. Emitters and detectors coupled at a distance of 3 to 4.5 cm allowed for the recording of 104 channels for each wavelength. Optodes were held in place perpendicularly, using a cap that was fitted on the head of participants in accordance with the 10 to 20 system.³¹ Optical intensity, including information regarding the average light intensity (DC), amplitude modulation (AC) and phase shift ($\varphi$), was measured with a sampling rate of 19.53 Hz (Boxy, ISS Inc., Champaign, Illinois). The channel setup covered both hemispheres equally and included the regions of interest for the investigation of language functions, i.e., frontal, temporal, and parietal lobes [Fig. 1(a)].

Fig. 1

fNIRS setup: (a) Probe placement (sources in small light gray dots, detectors in dark gray dots) as shown on an adult’s head model. (b) Experimental design including a 12-min resting-state followed by an expressive verbal fluency task that included 11 trials. Trials lasted 30 s and interstimulus intervals varied pseudo-randomly between 25 and 29 s.

Participants sat comfortably in a soundproof room. They were instructed to relax, to avoid any intentional movements or muscular tension, and to fix their gaze on the center of a screen placed at a distance of 114 cm. Participants underwent fNIRS recording during two conditions [Fig. 1(b)]: (1) a 12-min resting state with eyes open and (2) a verbal fluency task previously validated for expressive language-related activation.^32–34 The task consisted of 11 different familiar semantic categories (e.g., animals, colors, fruits, and so on) that appeared one at a time on the screen. Participants were instructed to name as many words as possible belonging to the specified category and to continue as long as the category name appeared on the screen. We used a block design paradigm in which periods of rest (fixation cross presented on the screen) and task (semantic category) alternated (Presentations^®, Neurobehavioral Systems, 2018). The interstimulus interval varied randomly between 25 and 29 s, while stimuli (one of the semantic categories) were always presented for 30 s. An audiovisual recording of the fNIRS session enabled visual support during offline preprocessing for the identification of movements and the timing of word articulation. Participants completed on average $10.4\pm 0.9$ blocks of the verbal fluency task and named an average of $14\pm 2.2$ words during each 30-s block. Data analysis was conducted with the use of a homemade toolbox (LIONirs)³⁵ adapted in SPM12 (Statistical Parametric Mapping)³⁶ in MATLAB^® (The MathWorks, Inc., Massachusetts).

2.2.

Validation Process

Validation of PARAFAC for artifact correction was done in two realistic applications as illustrated in Fig. 2. First, artifact correction with PARAFAC was tested on the real task-based signal of the whole sample, and its performance was compared with tPCA and ICA. Indices for the signal’s quality were used to compare the three techniques (Sec. 2.6). Second, a simulation analysis was conducted to investigate the efficacy of PARAFAC to correct artifacts with controlled parameters. Similarity metrics, quality measures and reconstruction of the ${\mathrm{HRF}}_{\text{sim}}$ were used to evaluate artifact correction with PARAFAC and tPCA both applied in a target manner (Secs. 2.6, 2.3, and 2.8).

Fig. 2

Processing streams to validate PARAFAC as a multidimensional artifact correction technique for fNIRS data. (a) The processing of the real task data derived from 17 subjects and (b) of the simulated data based on a single-subject, respectively. The metric(s) used to evaluate correction performance of each stream are presented on the right side. ${\mathrm{NIRS}}_{\mathrm{ini}}$ = initial signal used for comparison (with artifacts for the real task data; before artifact simulation for the resting-state data). PRD, Percent root difference; RMSE, root mean square error; $R$, Pearson product-moment correlation coefficient; SNR, signal-to-noise ratio; $R\lambda$, Pearson’s correlation between wavN elengths; GLM, general linear model; ${\mathrm{HRF}}_{\text{sim}}$, simulated hemodynamic response function.

2.3.

Nonsimulated Task-Related Motion Artifacts

Preprocessing of the task-condition data first included the automatic exclusion of channels with insufficient light intensity amplitude (average raw DC intensity across time $<100$). The signal was afterward segmented into blocks of 50 s (5 s resting-state baseline, 30 s task, and 15 s resting-state), and light intensity was converted to changes in optical density (normalization of each block). We then performed a semi-automatic artifact detection using a moving-window algorithm to automatically mark segments where an abrupt change of the signal’s variance exceeded three times the average variance of the previous interval³⁷ with a window duration of 0.8 s. Events that were 2 s or less apart were considered as one and channels that were strongly correlated with a noisy channel (Pearson correlation of $\ge 0.8$) within the aberrant segment were also marked as artifacted. The automatic detection step was reviewed afterward and adjusted based on an inter-rater visual inspection of the signal’s characteristics and the video recordings. Because the task required participants to name words aloud, artifacts in the current language paradigm were mostly due to facial movements related to the muscular contraction of articulation.

2.4.

Simulated Motion Artifacts

A real resting-state dataset from one of the participants in which we could identify a 180-s segment without any motion artifact was employed for the simulated artifact experiment. As for the task-based data, channels with insufficient raw light intensity amplitude were first excluded based on the previously mentioned criteria, yielding a total of 82 channels. The data were then converted to changes in optical density (normalization based on the whole 180-s segment). It served as the initial signal baseline (${\mathrm{NIRS}}_{\mathrm{ini}}$) for similarity and quality signal assessments.

To add a motion artifact with controlled parameters, we first used PCA to decompose all channels of a task-based fNIRS signal during a typical motion artifact (A) with a duration of 7 s. The obtained temporal signature of the first component was then added to the ${\mathrm{NIRS}}_{\mathrm{ini}}$. The spatial distribution of A, i.e., the weight of each channel was randomized, but the same distribution was added to both wavelengths with only a different overall scale. Figure 3 provides an overview of the set of simulations and their parameters, while more details are provided in Fig. S1 in the Supplemental Material. For the first simulation, the amplitude of the original artifact was modulated and scaled to five different amplitudes in order to produce artifacts with various signal-to-noise ratio (SNR) (simulations 1a to 1e). For the second simulation, we aimed to produce artifacts with more complex signatures (simulations 2a to 2f). The idea emerged from observations of artifact correction in the task-related signal and also aimed at exploring the decomposition when different artifacts are not completely orthogonal in the data. To do so, two individual artifacts (${\mathrm{A}}_1$, ${\mathrm{A}}_2$) with varying time intervals between the onset of ${\mathrm{A}}_1$ and ${\mathrm{A}}_2$ led to different temporal orthogonality ($r$) which was evaluated for each simulation. The levels of orthogonality between the time courses of the raw signal and the first artifact $r_a$ ($\text{Raw}\times A_1$) and between the time courses of the artifacted signal and the second artifact $r_b$ ($(\text{Raw}+A_1)\times A_2$), were derived from the angle between both signals (i.e., taking the time course of each signal as a vector of time points and computing the angle between both vectors). A normalized measure of the orthogonality level is then computed such that when the angle between the signals is 90 deg, the measure is maximum and equal to 1, while when the angle departs from 90 deg (both higher or lower) the measure linearly decreases to 0. Then, $r$ values of 1 indicate perfect orthogonality between the two time courses, while $r$ values lower than 1, specifically those closer to 0 indicate that the two time courses are not orthogonal to each other. Additionally, an HRF was synthesized (${\mathrm{HRF}}_{\text{sim}}$) in SPM by the linear combination of two gamma functions as proposed in the literature.³⁸ That is, gamma functions 1 and 2 had a time-to-peak of 5.4 and 10.8 s, respectively, a full-width-at-half maximum of 5.2 and 7.36 s respectively, and a scaling coefficient for the second gamma function of 0.35. The amplitude of the ${\mathrm{HRF}}_{\text{sim}}$ was scaled by 14% for the 830 nm signal and 0.6% for the 690 nm signal, which approximately produces an HRF with an $15\mu \mathrm{M}$ increase in HbO concentration and a $5\mu \mathrm{M}$ decrease in HbR concentration.^6,38–40 The response was simulated as a convolution of the HRF with a stimulus duration of 30 s resulting in a simulated physiological signal (${\mathrm{HRF}}_{\text{sim}}$) of $\sim 40\mathrm{s}$. The ${\mathrm{HRF}}_{\text{sim}}$ was then added to the signal along with an artifact ($A$), so as to produce different overlaps between both time courses (simulations 3a to 3d).

Fig. 3

Simulated fNIRS signals: (a) One example of the simulated fNIRS signal for each simulation condition one to three. (b) The details of the varying parameters of all simulations. Five scaling factors (0.4, 0.8, 1.0, 1.5, 3.0) were used to create simulated artifacts with different amplitudes. $\lambda 1|2=690|830\mathrm{nm}$; $A$, artifact; ${\mathrm{HRF}}_{\text{sim}}$ = simulated hemodynamic response function.

2.5.

Artifact Correction

In the task-based data set, artifact correction with PARAFAC and tPCA was applied in a target manner, i.e., to the segments identified during artifact detection as described in Sec. 2.3; while, with ICA the analysis was applied to the entire continuous signal. In the simulated signals, target correction was applied to a time interval of $\pm 2\mathrm{s}$ around the simulated artifacts. More details on each decomposition technique will be presented in the following sections (Secs. 2.5.1 and 2.5.2).

2.5.1.

Two-dimensional signal decomposition for artifact correction

For both ICA and tPCA, the decomposition of the 2D data matrix $X$ (whose elements $x_{i}j$ are indexed by channel $i$ and time points $j$), leads to $N_f$ components as defined in Eq. (1) and illustrated in Fig. 4(a)

Eq. (1)

$$x_{ij}=\sum _{f=1}^{Nf}{a}_{if}{b}_{jf}+{\mathcal{E}}_{ij}.$$

Fig. 4

Schematic representation of the decomposition models applied to fNIRS data ($A$). The data $X$ is arranged as a 2D data array by vertically concatenating the 2D matrices with dimensions being space ($s$) and time ($t$) for each wavelength ($w$). tPCA/ICA decomposes the array into components, each being a bilinear product of the loading vectors representing temporal ($a_{\mathrm{tf}}$) and spatial signatures ($b_{\mathrm{swf}}$). The latter is formed by the spatial signatures for the different wavelengths, which are represented in components without taking into account their spatial dependence, i.e., for the same temporal signature of each component, there will be two topographies corresponding to the two wavelengths. Matrices $A=\{a_f\}$ and $B=\{b_f\}$, contain as columns the temporal and spatial signatures for all components, respectively. ($B$) The data $X$ is arranged as a 3D data array with dimensions being time ($t$), space ($s$), and wavelengths ($w$). PARAFAC decomposes this array into the sum of components, each being a trilinear product of loading vectors representing temporal ($a_{\mathrm{tf}}$), spatial (channel, $b_{\mathrm{sf}}$), and spectral (wavelength, $c_{\mathrm{wf}}$) signatures. In practice, the decomposition consists of finding the matrices $A=\{a_f\}$, $B=\{b_f\}$, and $C=\{c_f\}$ that explain $X$ with minimal residual error.

The matrix concatenates two separate data matrices ($X_1,X_2$) corresponding to the 2D structure for both wavelengths. The maximum number of components $N_f$ could be equal to or less than the smaller dimension of the $X$ matrix dimension, i.e., twice the number of channels or number of time points. Each component is modeled as the product of two factors/vectors which represent signatures of the space ($a_f$) and time ($b_f$) dimensions. The temporal signatures are constrained to be orthogonal among components, and rotated in order to obtain those signatures that offer the highest variance explanation (varimax) from all the infinite solutions of the decomposition. The unexplained part of the data is considered irrelevant activity or noise ($\mathcal{E}$).

Artifact correction with tPCA was performed on each time interval containing artifacts as specified by the time interval of the simulated artifacts or the identified artifact detection of the task-based data set. From the obtained components we subtracted the first, which explains most of the variance based on the assumption that the highest variance in the data for each segment is assumed to be caused by the artifact.

Artifact correction with ICA, was done using BrainVision Analyzer (Brain Products GmbH, Gilching, Germany), and followed the steps proposed by Plank.⁴¹ First, optical density fNIRS data were exported and ICA was applied on the entire unsegmented signal, as for the method commonly reported in the literature.¹⁸ Noisy intervals reflecting artifacts were identified by semiautomatic inter-rater artifact detection (Sec. 2.3). We visually selected the ICA components reflecting the artifact, based first on their time course, i.e., high variation of light intensity similar to the artifact signature. Second, we rejected those whose subtraction would induce artificial artifacts elsewhere. Data were subsequently imported back into the LIONirs toolbox³⁵ in order to apply segmentation (Sec. 2.3).

2.6.

Signal Quality

The signal’s overall quality was estimated by two quality measures:^4,6,45

2.7.

Similarity Indices

For the analysis of simulated data, three additional metrics were used to evaluate performance of artifact correction: (1) the percent root difference (PRD), (2) the root mean square error (RMSE), and (3) the Pearson product-moment correlation ($R$). These allowed us to estimate the degree of correspondence between the initial signal without simulated artifacts (${\mathrm{NIRS}}_{\mathrm{ini}}$ and ${\mathrm{NIRS}}_{\mathrm{ini}}+{\mathrm{HRF}}_{\text{sim}}$) and the signal after correction of the simulated artifacts.^40,46,47 PRD, RMSE, and $R$ were defined as follows:

For all three indices, $x(t_i)$ and $y(t_i)$ represented the $i$’th point of the time courses of the corrected signal and of the initial artifact-free signal (${\mathrm{NIRS}}_{\mathrm{ini}}$ and ${\mathrm{NIRS}}_{\mathrm{ini}}+{\mathrm{HRF}}_{\text{sim}}$), respectively. $N$ corresponded to the duration of the time courses and $\overline{x}$ and $\overline{y}$ to their respective mean values along time. An ideal artifact correction would have been achieved when the artifact was completely removed, and the corrected signal maximally resembled the initial signal. As PDR and RMSE inform on the difference of the two signals, the smaller they were, the more accurate the performance of artifact correction was. On the contrary, the $R$ index represents the similarity of both signals, hence a higher value suggested better artifact correction.

2.8.

HRF Recovery

To investigate how artifact correction would affect the interpretation of the hemodynamic response, in the last simulations (3a to 3d) we aimed to recover a previously synthesized HRF. Therefore, a Butterworth low-pass filter (filter order = 4, cut-off frequency = 0.2 Hz) was applied to the signal to remove oscillations caused by heartbeat and respiration. The optical intensity changes of the two wavelengths were then transformed into relative concentration changes of HbO and HbR, done so by using an age-adapted differential pathlength factor (DPF) and the modified Beer–Lambert Law.^48,49 A general linear model (GLM) with the ${\mathrm{HRF}}_{\text{sim}}$ as a predictor was then applied. The GLM assumes that a linear relation exists between different inputs,⁵⁰ and has already been applied in several fNIRS studies to specify the HRF.^7,51,52 The cerebral activation of the signal ($y$) is defined as

2.9.

Statistical Analysis

Quality and similarity measures were computed for each channel separately and then averaged across channels for statistical analysis. Prior to statistical analysis, correlation coefficients ($R$ and $R\lambda$) were standardized using Fisher’s transformation to obtain values following a normal distribution. For the real task-based data, we conducted statistical analysis to compare the signal quality metrics among the motion artifact correction techniques. A repeated measures ANOVA including PARAFAC, tPCA, ICA, and the uncorrected reference signal as the within factor, was performed independently for the mean of each metric across all channels (i.e., SNR and $R\lambda$). In case of the $R\lambda$, the reference of the baseline (nonartifacted signal) was included as another condition. Follow-up paired contrasts were conducted with a critical alpha of 0.05. For the simulations, we applied a repeated measures ANOVAs to compare the mean outcome of the quality and similarity metrics between the different simulations, i.e., 1a to 1e for the amplitude scaling, 2a to 2f for the complex artifacts and 3a to 3d for the ${\mathrm{HRF}}_{\text{sim}}$, and within correction condition, i.e., PARAFAC, tPCA, and the uncorrected signal. Tukey correction for multiple comparisons was applied for posthoc analysis.

3.1.

Correction of Nonsimulated Task-related Motion Artifacts

First, decomposition with PARAFAC was applied for the correction of real task-based motion artifacts. The signal’s quality (SNR, $R\lambda$) after correction with PARAFAC was evaluated and compared to the signal before correction (uncorrected), the signal after correction with two 2D decomposition techniques (i.e., tPCA and ICA), and a reference signal of a segment without artifacts ($R\lambda$). Figure 5 provides an example of target artifact correction with PARAFAC illustrating the 3D decomposition for a motion artifact detected in the real task-based signal. The time course of the three components obtained with PARAFAC allowed to differentiate between the two components with distinct artifact signatures, and the one component representing the clean signal. Spatial distribution of each component is shown as a loading matrix and helps to identify in which channels (localization) the artifact’s signatures were most important. The two components that appear to represent the artifact showed similar scores for both wavelengths, as represented in the wavelength signatures. Subsequently, removal of PARAFAC first and second components resulted in a less noisy time interval and satisfyingly corrected signal.

Fig. 5

Example of target 3D PARAFAC decomposition to correct motion artifacts in a task-based fNIRS data set. (a) The initial uncorrected data segment. (c) The temporal, spatial (channel), and wavelength ($\lambda 1|2=690|830\mathrm{nm}$) signatures of the components identified with PARAFAC decomposition are presented, respectively. (c) The corrected signal, illustrating the efficacy of movement artifact correction with PARAFAC after subtraction of two components (1 and 2). $Y$-axis is presented in arbitrary units.

Automatic and manual artifact detection agreed on a majority of the signal with an average concordance of 92.8%, so only minor manual adjustments had been applied. A total of 585 different artifactual events, on average 66 events per subject with a mean, minimal and maximal duration of 5, 0.7 and 21.2 s, respectively, were considered for correction. Comparison of the SNR included measures of every channel of both wavelengths. Repeated-measures ANOVA with Greenhouse–Geisser correction for the SNR of all conditions (uncorrected, PARAFAC, tPCA, and ICA) revealed a significant main effect, $F(\mathrm{2.3,44966.3})=7550.7$, $p<0.001$, $\eta =0.28$. The uncorrected signal had the largest ratio of noise ($-0.69$), followed by the signal after ICA correction ($-0.50$), tPCA correction ($-0.33$), and PARAFAC correction ($-0.31$) (Table 1). Pairwise contrasts (Cohen, 1988) revealed that PARAFAC correction resulted in a significantly higher SNR compared with the uncorrected signal ($\eta =0.45$), and the signal after correction with ICA ($\eta =0.13$), both having a large effect size. Comparison of PARAFAC and tPCA, though statistically significant, revealed only a small effect size ($\eta =0.002$). Repeated-measures ANOVA with Greenhouse–Geisser correction for the $R\lambda$ revealed a significant difference between the five conditions: $F(\mathrm{3.2,30999.0})=2686.46$, $p<0.01$, $\eta =0.22$. Mean $R\lambda$ coefficients showed the highest association between wavelengths for the uncorrected signal (0.74), followed by the signal after ICA correction (0.67), tPCA correction (0.53), PARAFAC correction (0.51), and the artifact-free reference signal (0.47) (Table 1). Contrasts of $\mathrm{R}\lambda$ coefficients between PARAFAC and ICA ($\eta =0.14$, large effect) as well as the uncorrected signal ($\eta =0.35$, large effect) suggested that wavelengths were greatly less correlated after PARAFAC correction. Correction with PARAFAC resulted in slightly higher $R\lambda$ coefficients compared to segments without artifacts (artifact-free reference, $\eta =0.01$, small effect), and slightly lower $\mathrm{R}\lambda$ coefficients compared to correction with tPCA ($\eta =0.01$, small effect).

Table 1

Mean quality metrics for evaluation of artifact correction in the task-based data set.

3.2.

Correction of Simulated Motion Artifacts

As a second step, we evaluated the performance of PARAFAC to correct simulated artifacts with varying parameters, see Fig. 3 for the composition of the different scenarios. Three similarity indices (PRD, RMSE, and $R$) for the agreement of the initial (${\mathrm{NIRS}}_{\mathrm{ini}}$) and the corrected signal (${\mathrm{NIRS}}_{\mathrm{PARAFAC}/\mathrm{tPCA}}$), two quality measures (SNR, $R\lambda$) and the reconstruction of the synthesized HRF were used to quantitatively compare performance of PARAFAC to tPCA. The sample for each index (PRD, RMSE, $R$, SNR) consisted of 164 measures, corresponding to channels of both wavelengths. Since $\mathrm{R}\lambda$ is a correlation coefficient between both wavelengths, the sample consisted of 82 values. For each simulation condition (1a to 1e, 2a to 2f, and 3a to 3d), indices were computed for three correction conditions (=within-subject factor), i.e., the signal corrected with PARAFAC or tPCA and the uncorrected signal. Table 2 shows the results of the interaction effects between simulations and corrections as revealed by the repeated measures ANOVA. Results of post-hoc comparisons with Tukey correction are illustrated in Fig. 6 using the example of PRD. The results of the other indices are mostly in line with those. Detailed results for the RMSE, $R$, SNR, and $\mathrm{R}\lambda$ can be found in Figs. S2–S5 in the Supplemental Material, respectively.

Table 2

Interaction effects of the statistical comparisons between simulations and within correction methods.

Fig. 6

Evaluation of correction in simulated motion artifacts based on signal similarity. Simulations are identified on the $x$-axis. [(a) 1a–1e] artifacts with different amplitude sizes; [(b) 2a–2f] complex artifacts with two superimposed artifacts and an onset delay between the first ($A_1$) and second artifact ($A_2$); [(c) 3a–3d] the onset of the artifact relative to the beginning of a simulated HRF. Performance of artifact correction is illustrated by the use of the similarity index percent root difference (PRD), where a lower value represents higher resemblance between the corrected and the initial clean fNIRS signal, hence a better correction of the artifact. Results are displayed separately for the correction with PARAFAC (light gray bars) and tPCA (gray bars). The PRD of the uncorrected signal is not displayed in this figure but differed significantly from both correction techniques in all conditions, except where specified (n.s.) otherwise inside the bar. Significance level are based on post-hoc tests with Tukey correction. ***$p\le 0.001$, n.s. $p>0.05$. $\text{Uncorrected}={\mathrm{NIRS}}_{\mathrm{ini}}$ + artifact ($A_1$) without correction, $\mathrm{PARAFAC}={\mathrm{NIRS}}_{\mathrm{ini}}$ + artifact ($A_1/A_2$) after artifact correction with PARAFAC, $\mathrm{tPCA}={\mathrm{NIRS}}_{\mathrm{ini}}$ + artifact ($A_1/A_2$) after artifact correction with tPCA.

Simulations 1a to 1e included artifacts with extra small, small, medium, large, and extra-large amplitudes. Similarity indices (PRD, RMSE, and $R$) revealed that both PARAFAC and tPCA resulted in a significant higher resemblance of the corrected signal and the initial signal compared to the uncorrected condition for all amplitude sizes except for the artifact with an extra small amplitude, where tPCA did not lead to a significant improvement compared with the uncorrected signal. Similarly, both correction methods led to a significantly higher signal quality, i.e., lower SNR and $R\lambda$, compared with the uncorrected signal with the same exception for tPCA in the extra small amplitude, where no difference was observed as compared to the uncorrected signal. When artifacts had a small, medium or extra-large amplitude, correction methods obtained comparable results over all indices and did not statistically differ. The most consistent differences between correction methods can be observed for the artifact with extra small and large amplitudes. While similarity metrics and quality measures revealed that PARAFAC achieved significantly better results than tPCA for the correction of artifacts with an extra small amplitude, for artifacts with a large amplitude tPCA seemed to have a slight but statistically significant advantage as compared with PARAFAC.

Simulations 2a to 2f led to signals with different onset delays (1–5 s and 10 s) between two superimposed artifacts, which allowed to create different complex artifacts. Over all conditions, notwithstanding their onset delay, similarity indices indicated a significantly higher overlap with the initial signal for the signal corrected with PARAFAC as compared with the one corrected with tPCA and the uncorrected signal. Similarly, quality metrics revealed a significantly higher signal quality after correction with PARAFAC compared to the correction with tPCA and the uncorrected signal. Nevertheless, tPCA also resulted in statistically significant better results as compared to the uncorrected signal in all conditions both with regard to the similarity metrics and the quality measures.

In simulations 3a to 3d the onset of an artifact varied relative to the beginning of a simulated ${\mathrm{HRF}}_{\text{sim}}$: $+0$, $+5$, $+15$, and $+30\mathrm{s}$. Independent of the onset of the artifact, correction with PARAFAC resulted in a signal that showed a significantly higher overlap with the initial signal compared with the uncorrected condition. Similarly, quality measures showed a significant improvement of the signal’s quality after correction with PARAFAC compared with the uncorrected signal. tPCA achieved almost identical results as PARAFAC in comparison to the uncorrected signal, except when the artifact was placed shortly after the beginning of the ${\mathrm{HRF}}_{\text{sim}}$ ($+5\mathrm{s}$) for which PRD and RMSE did not show a significant difference between the uncorrected signal and the one corrected with tPCA. Further, for artifacts at the very beginning (+0 s), shortly after ($+5\mathrm{s}$) or at the very end ($+30\mathrm{s}$) of the ${\mathrm{HRF}}_{\text{sim}}$ (at the same time than the low-frequency increase or decrease that is intrinsic to the ${\mathrm{HRF}}_{\text{sim}}$), similarity indices revealed a significant better overlap of the corrected signal with the initial signal after correction with PARAFAC as opposed to tPCA. In contrast, tPCA compared to PARAFAC seemed to have a significant advantage for the correction of artifacts during the plateau of the ${\mathrm{HRF}}_{\text{sim}}$ ($+15\mathrm{s}$), which might be explained by a higher orthogonality between the artifact and the signal at this onset. Quality measures however only partially support these findings, namely in favor of correction with tPCA for the artifact at the beginning ($+0\mathrm{s}$) and in favor of PARAFAC for the artifact shortly after ($+5\mathrm{s}$) the beginning of the ${\mathrm{HRF}}_{\text{sim}}$. When the artifact is at the very beginning of the ${\mathrm{HRF}}_{\mathrm{sim}}$ ($+0\mathrm{s}$), signal quality is however significantly higher after correction with tPCA as compared with PARAFAC. There is no statistically significant difference between signal quality for any of the correction methods when the artifact is placed at the very end of the ${\mathrm{HRF}}_{\text{sim}}$ ($+30\mathrm{s}$). Further, results of the GLM revealed that in none of the conditions the simulated artifact led to an important reduction of the variance explained by the ${\mathrm{HRF}}_{\text{sim}}$ as compared to the initial signal (Table 3). In line with the results of the similarity metrics, reconstruction of the HRF after correction with PARAFAC was qualitatively better as compared with tPCA when the artifact occurred at the beginning ($+0\mathrm{s}$), shortly after ($+5\mathrm{s}$) and at the end ($+30\mathrm{s}$) of the HRF. In particular, for the artifact 5 s after the onset of the HRF, reconstruction after correction with tPCA seemed disturbed. There was no important difference between correction methods for the reconstruction of the HRF when the artifact was during the plateau of the HRF ($+15\mathrm{s}$).

Table 3

Percentage of explained variance by the HRFsim based on the R2 of the GLM.

Orthogonality measure for the relation of the time courses of the raw signal and the artifact signal ($\text{Raw}\times A_1$: $r_a$) were calculated for all simulations, orthogonality for the time courses of the raw signal with the first artifact and the second artifact ($\text{Raw}{A}_1\times A_2$: $r_b$) was only computed for simulations of complex artifacts (2a to 2f). A value of 1, represented maximum orthogonality between both signals meaning they are in a 90 deg angle to each other. Any value smaller than 1 or even 0 indicated nonorthogonality meaning that they are in any other angle than 90 deg, maximally 0 deg or 180 deg to each other. Results of $r_a$ revealed consistently high orthogonality for all simulations with varying amplitudes ($\mathrm{Min}=\mathrm{Max}=0.97$) and with complex artifacts ($\mathrm{Min}=0.97$, $\mathrm{Max}=0.98$). Orthogonality for the simulations with the HRF revealed more variations ($\mathrm{Min}=0.88$, $\mathrm{Max}=0.97$) where the artifact 30 s after the HRFsim created the least orthogonal condition (Median ± standard error: $0.88\pm 0.03$), followed by the simulation with the artifact at 0 s (Median ± standard error: $0.93\pm 0.05$), 5 s (Median ± standard error: $0.96\pm 0.01$) and 15 s (Median ± standard error: $0.97\pm 0.01$) of the ${\mathrm{HRF}}_{\text{sim}}$. Results of the second orthogonality measure $r_b$ indicated that the onset delay of 2 s led to the least orthogonal condition (median ± standard error: $0.75\pm 0.07$), followed by the delay of 3 s (median ± standard error: $0.93\pm 0.02$) and 5 s (median ± standard error: $0.98\pm 0.02$). When the second artifact had a delay of 1, 4, or 10 s, there was equally high orthogonality for all three conditions (Median ± standard error: $0.99\pm 0.02$).

4.1.

Usefulness and Limitations of PARAFAC

PARAFAC is a data-driven approach based on the linear relations of the three dimensions of the fNIRS signal.^22,27,28 PARAFAC has the advantage of allowing for the three-way arrays of data to be uniquely decomposed into a sum of components, each of which is a trilinear combination of factors or signatures. The only statistical requirements of PARAFAC is that of a moderate linear independence across components, i.e. their time course, topography and wavelength characteristics. This is a less stringent requirement than previous models that underlie space/time decompositions (PCA or ICA). Each component provides characteristics of a particular pattern identified within the mixed measured fNIRS signal. When there are empirical or theoretical reasons to expect more than one relevant component (e.g., resting state functional connectivity analysis, epileptic activity, physiological aspects, and so on), PARAFAC would also allow disentangling several components. Moreover, when there are reasons to include other constraints such as orthogonality, non-negativity, smoothness and sparseness of the signatures, it has been shown that PARAFAC can also include them in the decomposition procedure.^22,57,58

Different from ICA, PARAFAC components can be ordered according to their importance in explaining the data’s variance.⁵⁹ This is why many studies using ICA apply data reduction or clustering techniques such as PCA prior to the execution of ICA.^18,60 Moreover, when using PARAFAC, the description of the data is based on more dimensions, and with less theoretical constraints (orthogonality and independence), which lead to the identification of fewer relevant components.^13,28,61,62 Compared with ICA, the selection of the relevant components is thus simplified. Importantly, when using PARAFAC, it is not recommended to ask for more than five components, as this increases the risk of overfitting the data with the decomposition model. Though the issue of selecting the optimal number of components should be better explored and standardized in future studies, based on our experience and previous studies, usually, three to five components sufficiently describe the signal’s relevant characteristics.^{13,22,23,28,29} The Corcondia index provides further information regarding a good fit of the PARAFAC model and the analyzed data. Precisely, a satisfying decomposition should have a core consistency of 85% or higher and usually drops below 80% if an extra component is extracted, suggesting overfitting.⁴⁴ Similar to PCA, a typical choice made when using PARAFAC analysis, is that of ordering the extracted components according to their contribution to explaining the variance of the data. In this sense, it can be expected that the first PARAFAC component will always represent the highest activity, which, in the case of target artifact correction, will correspond to the artifact activity. This makes PARAFAC also a promising choice to explore the development of simple automatic methodologies for the detection and correction of such artifacts. Relevant components for artifact correction could for instance be selected based on the amplitude size or the differences between scores for the two wavelengths. Despite being a time efficient and promising approach, we believe that at the current state visual verification cannot be ceased, because depending on the paradigm, noise may have characteristics that render automatic processing inadequate. Future studies should be devoted to this specific aim.

Some challenges of PARAFAC have already been discussed in previous studies.^13,28 One limitation is related to PARAFAC’s assumption of linear relations of the temporal characteristics between different channels [which also applies to (t)PCA and ICA]. Despite being one of the most common and simple models, imposing linearity may not entirely reflect complex cerebral processes.⁶³ Some attempts have been made to introduce approaches that tolerate nonlinear relations.^64,65 However, to date, linear models remain the most popular and adequate approach for analysis of macroscopic data in neuroscience. Thus, it represents a general limit of the domain, rather than one specific to PARAFAC. Another potential limitation of our implementation of PARAFAC might arise from the specific use of the simple alternating least-squares algorithm to perform the decomposition. Other methods have been proposed for estimating PARAFAC models, but some studies have shown that, given the uniqueness of the solution, they do not usually outperform the simple least-squares technique.⁶⁶

The need for careful preprocessing of data before characterizing brain activity constitutes another limitation. When PARAFAC was introduced for the analysis of EEG data, the authors emphasized the importance of searching for constant factors, outliers, and degeneracy.^13,26 Detailed information on how to deal with these aspects can be found in previous literature.¹³ In the current study, the implementation of PARAFAC into the LIONirs toolbox,³⁵ gave us a certain flexibility regarding the channels and time segments to be included. Thus, we were able to exclude deviant channels prior to use PARAFAC.

Finally, among the applications of PARAFAC not yet tested in fNIRS, is its use to identify patterns of brain activity in the hemodynamic response as it has been applied in EEG data analysis.^13,29,42 Compared with other techniques such as GLM or global averaging who treat each channel independently and provide a narrower spatial representation of the dominant activation, PARAFAC analysis could reveal a wider distributed activation pattern. This would strongly correspond to the current understanding of cerebral processing, whereby mostly large-scale networks, and not isolated regions, are considered to be involved in various cognitive processes.^55,67,68 PARAFAC could thus appear to be suitable for reflecting cerebral processes occurring in distributed networks, rather than for the identification of a specific core region. This is obviously given a sufficient fNIRS covering. Moreover, the topographic signature would correspond to a whole time course of the extracted activations and can easily be subjected to diffusion optical tomography in order to locate the HbO and HbR concentration changes in the brain cortex.⁶⁹ Another potential application is using PARAFAC as a screening tool. For instance, Miwakeichi and colleagues¹³ applied PARAFAC to an EEG dataset in order to extract one component related to ocular movement artifacts, and subsequently used the PARAFAC analysis fixing spatial and spectral signatures of that component to screen a second dataset in order to identify and correct similar artifacts. This application can also support the development of a detection method as has already been done in other fields such as detection of epileptic seizures^29,70 and for brain–computer interfacing.^71–74 Although this can be useful for correcting artifacts that have consistent topographical and wavelength profiles, movement artifacts were mostly related to articulation in our study, often showing quite different signatures. A different paradigm may enable the testing of this application of PARAFAC in fNIRS.

Our findings also encourage the general application of PARAFAC for multidimensional decomposition in neuroimaging, where data can often be described in more than two dimensions. PARAFAC could for instance be of interest for fNIRS data acquired with a multiwavelength ($>2$ wavelengths) system.^75,76 Beyond the obvious dimensions of a technique such as that of time × space × wavelength (fNIRS) or frequency (EEG), characteristics of the task paradigm or group variables could also be considered in the decomposition model.⁴² For instance, if a “group” dimension was to be included, groups could be compared regarding their relative weight over the main cerebral activation component that would have been outlined from the PARAFAC model.

4.2.

Limits of the Current Study

The use of a language paradigm allowed for the induction of artifacts that were mainly related to articulation and often overlapped with the stimulus onset. This makes identification and correction of artifacts particularly challenging and may limit the generalization of our results. From video recordings, however, we were able to acknowledge that swallowing, jaw or tongue movements, eye blinks, and frowning also induced artifacts in our data set. Our results particularly contribute to an improved artifact correction in a paradigm where artifacts can pose a significant problem for interpretation. The validation of PARAFAC in a realistic scenario represents an important advantage because simulated artifacts either completely artificial or by instruction of participants to make certain movements, rarely live up to the complexity and variability of real spontaneous artifacts.⁵ Though the simulated artifacts were also based on one identified during the task-condition, parameters were modified and allowed to create noise with different characteristics. Further fNIRS studies using PARAFAC, and including different task paradigms, resting-state, and recording in naturalistic environments when participants are moving or are instructed to perform specific movements, could show its utility for a wider range of artifacts (e.g., baseline shifts or physiology).

The validation of PARAFAC for artifact correction was conducted by comparing results with two other decomposition techniques. This could be considered one limit of the present study, because we did not consider other correction approaches for comparison. The objective of this study was not to conduct exhaustive and systematic comparisons in order to identify the best method for artifact correction, but rather to introduce PARAFAC as an alternative and adequate method in fNIRS analysis. Future studies could therefore expand the validation of PARAFAC and extend comparisons of its performance to other tools that do not necessarily use decomposition, such as spline interpolation or wavelet filtering. With relation to that, even if there were to be a slight advantage of PARAFAC over tPCA in artifact correction results, the small difference might not have a concrete impact on the applicability of those two methods. Qualitative differences between tPCA and PARAFAC regarding the type of artifacts and corrected signal should be further assessed and described in subsequent studies to better understand how their efficacy may vary according to the artifact’s characteristics.

Last but not least, PARAFAC’s components related to artifacts were selected during a visual inspection based on their temporal overlap with the artifact, i.e., sudden change of amplitude, and those showing similar weights for both wavelengths. The smallest number of possible components that sufficiently corrected the artifact were selected. Even though our results show that this procedure led to satisfying correction, it is not a fully standardized procedure, which may hamper an easier, automatic correction that would reduce the need of the researcher’s expertise in detecting and adequately applying correction. This issue could be addressed using the screening procedure mentioned in the previous section or using other templates for the artifacts of interest or by applying the usual automatic rules based on amplitude characteristics to the temporal signatures obtained by PARAFAC instead of the original mixed data. The goal of the current study was not to propose a method for artifact detection but to validate the relative efficacy of these methods in artifact correction. It would be interesting for future studies to address the development of a more automatic procedure based on PARAFAC analysis contributing to a faster and standardized processing pipeline.