A multi-step blind source separation approach for the attenuation of artifacts in mobile high-density electroencephalography data

Electroencephalography (EEG) is a technique that measures electrical potentials over the scalp, which vary over time as a result of neuronal activity in the brain. EEG signals, however, are not only related to neural sources, but also reflect several kinds of other sources: blinks, heart-beat, muscle contractions, movements, as well as environmental noise (Urigüen and Garcia-Zapirain 2015). It is important to note that the impact of non-neural signals in EEG data, also referred to as artifacts, can be particularly strong in the presence of movements and muscle contractions of the participant. This currently limits the use of the EEG technology in naturalistic environments. Accordingly, EEG is often used for studies that can be conducted under very controlled experimental conditions (Dietrich and Kanso 2010, Liu et al 2017, 2018, Samogin et al 2019, 2020, Zhao et al 2019).

A considerable number of methodological studies have focused on automated EEG artifact attenuation. Traditionally, artifacts are estimated and linearly regressed from EEG signals (Gratton et al 1983, Kenemans et al 1991, McMenamin et al 2009), or are adaptively filtered (Narasimhan and Dutt 1996, He et al 2004, Correa et al 2007, Jafarifarmand and Badamchizadeh 2013). Those methods require the availability of external signals from which artifactual contribution can be estimated, such as electrooculography (EOG) and electromyography (EMG) signals. The artifacts can also be attenuated without the use of external signals, using methods relying on wavelet transformation. Those methods first decompose the EEG data using wavelet coefficients, and then reconstruct new signals having decreased artifactual coefficients (Zikov et al 2002, Krishnaveni et al 2006, Bajaj et al 2020). Their performance varies with the choice of mother wavelets, and the criteria used for artifactual coefficient attenuation. Another group of methods uses empirical mode decomposition (EMD) to decompose the EEG data into narrow-band intrinsic mode functions (IMFs), and to reconstruct signals by mixing the IMFs that are not deemed to be artifactual (Zeng et al 2013, Chen et al 2014, Sweeney-Reed et al 2018). Notably, the performance of EMD approaches may be limited by the mode mixing problem (Zheng et al 2014).

In recent years, the development of the high-density EEG (hdEEG) technology has permitted to attain a finer spatial sampling of scalp potentials (Michel and Murray 2012). This allowed a reliable use of blind source separation (BSS) for EEG artifact attenuation (Delorme et al 2012, Hu et al 2018). BSS methods linearly decompose input signals into a set of components that are linearly mixed over the channels. Such components can be retrieved using different BSS methods. A classical BSS approach is principal component analysis (PCA), which yields components that are orthogonal (i.e. uncorrelated) to one another (Wold et al 1987). Canonical correlation analysis (CCA) is a variation of the PCA, which enforces the components to be maximally autocorrelated and mutually uncorrelated (De Clercq et al 2006). Furthermore, independent component analysis (ICA) approaches decompose input signals into statistically independent components (Hyvärinen and Oja 2000, Jung et al 2000). Different ICA implementations have been developed over the years, among which the most widespread are: fast fixed-point independent component analysis (FastICA) (Hyvarinen 1999), Information maximization (Infomax) (Makeig et al 1996), and joint approximate diagonalization of eigenmatrices (JADE) (Cardoso and Souloumiac 1993). More recently, independent vector analysis (IVA) has been proposed as an extension of ICA, in which the criterion of component independence is preserved, but applied to multiple datasets derived from the input data; at the same time, maximum dependence is required between components extracted from different datasets (Anderson et al 2011, Chen et al 2017). Generally speaking, when BSS methods are applied to EEG data, the output components are classified into two categories: artifactual and non-artifactual (Liu et al 2017). Then, neural signals are obtained by re-mixing the non-artifactual components using the associated weights, also produced by BSS (Liu et al 2017).

A key point of using BSS for EEG artifact removal is the general superiority in the reconstruction of artifactual components, as compared to alternative approaches (Pion-Tonachini et al 2019). Several studies using a manual classification of BSS components performed through visual examination of their temporal, spatial and spectral features (Chaumon et al 2015). However, the classification can also be performed in an automated or semi-automated manner, using tools such as FASTER (Nolan et al 2010), ADJUST (Mognon et al 2011) and ICLabel (Pion-Tonachini et al 2019). These tools have been developed exploiting ICA but can be in principle used with any BSS method. Notably, it has been suggested that different BSS methods are suited to isolate specific kinds of artifacts. For instance, ICA has been consistently found to be effective when EEG data are corrupted by ocular artifacts (Jung et al 1998, 2000); conversely the use of IVA has been found to be superior to ICA and CCA for the attenuation of myogenic artifacts (Chen et al 2017). To the best of our knowledge, previous studies always applied a single BSS approach for the attenuation of artifacts mixed in EEG data (Nolan et al 2010, Mognon et al 2011, Pion-Tonachini et al 2019). Furthermore, most studies focused on the ocular artifact, as this is the most prominent in highly controlled EEG experiments (Kenemans et al 1991, He et al 2004, Hsu et al 2012, Urigüen and Garcia-Zapirain 2015, Guarnieri et al 2018). When several kinds of artifacts are mixed in the EEG data, as for instance in mobile brain/body imaging studies (Nathan and Contreras-Vidal 2016, Jungnickel et al 2019), the use of a single BSS method is likely to be suboptimal (Barban et al 2021).

In this study, we propose a novel multi-step BSS approach for attenuating ocular, movement and myogenic artifacts, which are typically present in mobile hdEEG data (Urigüen and Garcia-Zapirain 2015, Nathan and Contreras-Vidal 2016, Richer et al 2020). To optimize and validate our novel approach, we first conducted an event-related potential (ERP) analysis (Okita 1979) of hdEEG data collected in either standing, slow-walking or fast-walking conditions, during auditory stimulation. Then, we used an additional hdEEG dataset collected in participants who alternated either slow-walking or fast-walking periods with stance periods, to study movement-related modulations of neural oscillations by performing an event-related desynchronization/synchronization (ERD/ERS) (Pfurtscheller and Da Silva 1999) analysis. ICLabel (Pion-Tonachini et al 2019), a single-step BSS approach widely used for artifact removal, was used for benchmarking purposes. We also included in the performance analysis the single-step BSS method previously used by our group (Liu et al 2017). Specifically, we tested whether the multi-step BSS approach, by more effectively attenuating the artifacts mixed in the data, can yield stronger and more reliable neural responses than alternative solutions.

2.1. Description of the multi-step BSS approach

Any BSS method models a set of observations $X\left( t \right)$ as a linear combination of source signals $S\left( t \right)$, represented through the mixing matrix A (Cao and Liu 1996). The relationship between source signals and observations is expressed by the following formula:

$$\begin{equation}X\left( t \right) = A \cdot S\left( t \right).\end{equation} \tag{ 1 }$$

Starting from the observations, each BSS method produces estimates of the source signals $\hat S\left( t \right)$ decomposed using different mathematical criteria with specific numerical implementation. In general, source signal estimates can be expressed as the linear combination of the observations, represented through the unmixing matrix W:

$$\begin{equation}\hat S\left( t \right) = W \cdot X\left( t \right).\end{equation} \tag{ 2 }$$

In order to use BSS for EEG artifact removal, the source signals need to be classified in two categories: neural or non-neuronal. Notably, the accuracy of source signal classification is dependent on the reliability of source signal reconstruction by BSS (Fitzgibbon et al 2007). In turn, this depends on the typology and magnitude of the artifacts that are mixed in the EEG data (Fitzgibbon et al 2007). Once the sources are classified, it is possible to create a spatial filtering matrix $Z$. This is a diagonal matrix, with the diagonal element ${Z_{ii}}$ being equal to 1 if the source is classified as neuronal, and equal to 0 otherwise (Mantini et al 2008). The artifact-corrected EEG data can then be obtained starting from the original EEG data by using the following formula:

$$\begin{equation}{X_c}\left( t \right) = A \cdot Z \cdot W \cdot X\left( t \right).\end{equation} \tag{ 3 }$$

The formula above is classically used when a single BSS method is applied to the EEG data for artifact removal (Mantini et al 2007a, 2008). In this case, ICA is very often used as BSS method, with the Infomax algorithm (Bell and Sejnowski 1995) and the FastICA algorithm with either deflation or symmetric decomposition (FastICA-defl/symm) (Hyvarinen 1999) being the most widespread solutions.

We mainly focused on the development of a method for the attenuation of ocular, movement and myogenic artifacts, which are the typical artifacts that can be found in mobile hdEEG data (Urigüen and Garcia-Zapirain 2015, Nathan and Contreras-Vidal 2016, Richer et al 2020). Specifically, we devised a three-step BSS solution to remove those artifacts in a sequential order (figure 1). We first removed ocular artifacts, then movement artifacts and finally myogenic artifacts, as we gave precedence to the attenuation of artifacts with higher relative power (figure S1, available online at stacks.iop.org/JNE/18/066041/mmedia) (Urigüen and Garcia-Zapirain 2015, Nathan and Contreras-Vidal 2016, Richer et al 2020). The relative power has been used as criterion for defining the order of the artifact removal steps, since BSS methods are sensitive to the signal-to-noise ratio (SNR) in the input data (Fitzgibbon et al 2007). For each artifact removal step, the attenuation of artifacts with highest relative power increases the SNR of the residual artifacts, and hence their separation by BSS in the subsequent steps.

Zoom In Zoom Out Reset image size

For each step, we optimized our approach by selecting the BSS method that provided the best decomposition performance, among the following seven methods: PCA (Wold et al 1987), CCA (De Clercq et al 2006), FastICA-defl/symm (Hyvarinen 1999), Infomax (Makeig et al 1996), JADE (Cardoso and Souloumiac 1993) and IVA-G (Anderson et al 2011). The performance was defined using the maximum absolute correlation between estimated source signals and reference artifactual signals (e.g. those derived from EOG, EMG and kinematic signals).

For each step, we also defined a feature to classify the artifactual sources. The ocular artifacts were detected using the maximum of the signal kurtosis calculated over contiguous time windows (Balanda and MacGillivray 1988); the movement artifacts were isolated by means of the mean of the sample entropy calculated over contiguous time windows (Richman et al 2004); the myogenic artifacts were classified using the ratio between the signal power at high frequencies as compared to the total power (Urigüen and Garcia-Zapirain 2015). We optimized the window size (for ocular and movement artifacts) or the cut-off frequency (for myogenic artifacts), as well as the specific classification threshold using a receiver operating characteristic (ROC) analysis (Mantini et al 2008). To this end, the source signals obtained using BSS were manually and independently classified by two experts. The source signals were then split into a training and a testing dataset, each including six participants. Using the training dataset, we selected the parameter configuration (window size/cut-off frequency) that maximized the area under the ROC curve (AUC), and correspondingly defined the threshold for which type I and type II errors were balanced (Mantini et al 2008). The testing dataset was then used to quantify, in an independent manner, the classification accuracy, defined as the relative number of correctly classified sources.

2.2. Validation of the proposed approach

2.2.1. Experiment

Twelve right-handed participants (seven females and five males, age 22–31 years) were recruited for this study. All participants were not affected by any brain-related disease or other medical condition. The experiment, which was approved by the ethical committee of the Liguria Region (reference: 238/2019), was conducted in accordance with the Helsinki declaration and its later amendments. An informed consent was acquired from each participant before the experiment.

The experiment was divided in two parts. The first one comprised three conditions lasting 6 min in total: standing, slow-walking and fast-walking, whose order was counterbalanced across participants. In the slow-walking and fast-walking conditions, the participants turned around after a distance of about 8 m. Periods of 2 s preceding and following turns, respectively, were removed from the analyses. The walking speed was chosen by each participant and resulted in 1.7–2.2 km h⁻¹ and 2.8–3.3 km h⁻¹, for slow and fast walk, respectively (figure S2(a)). During each condition, a series of auditory stimuli were played. Each one consisted of a 900 Hz sinusoid lasting for 62 ms in total, with fade-in and fade-out periods of 10 ms and loudness of about 65 dB at the participant's ears. The inter stimulus interval was set to randomly vary between 1000 and 1375 ms (Zink et al 2016). Each condition included 320 auditory trials in total. Subsequently, in the second part of the experiment, the participants were asked to repeat the slow- and fast-walking conditions, but alternating periods of stance and movement of 6 s each, for a total number of 30 blocks (Zhao et al 2019). The walking speed was 1.6–2.0 km h⁻¹ and 2.6–3.5 km h⁻¹, for slow and fast walk, respectively (figure S2(b)). In this case, no auditory stimulation was delivered. The two conditions were counterbalanced across participants. Also in this case, the participants turned around after a distance of about 8 m, and periods of 2 s preceding and following turns, respectively, were removed from the analyses.

2.2.2. Data collection and preprocessing

During the experiment, we collected 128-channel hdEEG signals, 2-channel EOG signals and 2-channel EMG signals at 1 kHz sampling rate, using an ActiCHamp amplifier (Brain Products GmbH, Germany). The hdEEG sensors were embedded in an EEG cap and connected to the participant's scalp using electrode gel. All the EEG channels were referred to a single reference channel (identified as FCz, according to a standard 10/20 montage). Bipolar EOG sensors were used to record horizontal and vertical eye movements. Bipolar EMG sensors were placed bilaterally over the neck to record muscular activity of left and right splenius capitis. The amplifier was positioned in a backpack worn by the participant. An auxiliary amplifier channel was connected to a microphone, for recovering precise sound stimulus timing. In parallel to hdEEG/EOG signals, we also collected acceleration data at 148 Hz sampling rate using a Trigno wireless system (Delsys Inc., USA). A total of 14 acceleration sensors were bilaterally positioned over the participant's limbs, in correspondence of the following body parts: arm, wrist, anterior thigh, anterior leg, posterior thigh, posterior leg, and ankle (Artoni et al 2017). The signal from the microphone was also sent to a Trigno wireless sensor, for synchronization purposes. Just after the experiment, we obtained a 3D scan of the participant's head using an iPad (Apple Inc. Cupertino, CA, US) equipped with a Structure Sensor (Occipital Inc., Boulder, CO, US).

In a separate session, the structural MR image of each participant's head was collected with a 3 T Philips Achieva MR scanner (Philips Medical Systems, the Netherlands) using a T1-weighted magnetization-prepared rapid-acquisition gradient-echo sequence. The scanning parameters were: repetition time (TR) = 9.6 ms, echo time (TE) = 4.6 ms, 160 coronal slices, 250 × 250 matrix, voxel size 0.98 × 0.98 × 1.2 mm³.

EEG electrode locations of each participant were obtained from the 3D scan of the head using the SPOT3D software (Taberna et al 2019a, 2019b). HdEEG, EOG, EMG and acceleration data were digitally filtered in the band 1–80 Hz, in line with our previous studies (Mantini et al 2007b, Liu et al 2017, Samogin et al 2019, 2020, Zhao et al 2019). Furthermore, hdEEG data were corrected for bad channels (Guarnieri et al 2021) and re-referenced to average reference (Liu et al 2015).

2.2.3. Performance analysis

Artifacts in mobile hdEEG data were attenuated using our novel multi-step BSS approach, as well as two validated one-step BSS methods, which were used as benchmark. The first of the one-step BSS methods was developed by our group, and entailed data decomposition by FastICA-defl and artifact detection by using a threshold-based classification approach (Liu et al 2017). The second one-step BSS method was based on the Infomax algorithm for data decomposition and on ICLabel, a multivariate classifier, for artifact detection (Pion-Tonachini et al 2019).

We first conducted in-depth analyses concerning the multi-step BSS approach. Specifically, we evaluated the EEG traces before and after each artifact attenuation step, for each condition. We quantified the reduction in artifactual contribution using maximum absolute correlations between hdEEG and EOG signals for the first step, between hdEEG and movement velocity signals for the second step, and between hdEEG envelopes and neck EMG envelopes for the third step. Envelopes were obtained as the absolute value of the Hilbert transform of the signals.

We also extracted performance measures not only from data processed with multi-step BSS, but also unprocessed and processed with one-step BSS (either FastICA-defl or Infomax). These measures entailed the ERP for the auditory part of the experiment (Zink et al 2016), and the ERD/ERS for the motor part of the experiment (Pfurtscheller and Da Silva 1999, Zhao et al 2019).

ERP time courses were obtained by epoching the hdEEG data with a [−100, 400] ms window with respect to auditory stimulus onset, and then averaging the resulting epochs. For each ERP time-course, we estimated the noise level as the averaged standard deviation in a [−100, 0] ms window, and then averaged the resulting noise values across all channels. In line with previous auditory ERP studies, the signal level was quantified as the mean absolute intensity in a [80, 120] ms window for the ERP time-course at channel Cz (Oray et al 2002, Miller and Zhang 2014). We also assessed the topographies of the auditory N1 response, for latencies ranging between 80 and 120 ms. Grand-averages for ERP time courses and topographies were obtained by averaging data across participants.

ERD/ERS time-courses and maps were obtained from hdEEG data collected during slow- and fast-walking (figure S2(a)), using the procedure described in one of our recent studies (Zhao et al 2019). To this end, we performed hdEEG source localization through the following three steps. (1) Head tissue segmentation. We segmented the individual MR images into 12 tissue layers (skin, eyes, muscle, fat, spongy bone, compact bone, cortical/subcortical gray matter, cerebellar gray matter, cortical/subcortical white matter, cerebellar white matter, cerebrospinal fluid, and brain stem) using the MR-TIM software (Taberna et al 2021a). The conductivity of each tissue was set according to previous studies (Haueisen et al 1997, Holdefer et al 2006). (2) Leadfield matrix calculation. The leadfield matrix, which expresses the relationship between neural activity and electrical potentials over the scalp, was calculated using the Simbio software (Vorwerk et al 2018), integrated in the FieldTrip toolbox (Oostenveld et al 2011). In detail, we generated hexahedral meshes from the individual 12-tissue layer MR images, and the elements of the meshes in the gray matter were set to represent possible sources of neural activity. The leadfield matrix was then interpolated to a 6 mm isotropic volumetric grid. (3) Neural activity reconstruction. We fed the leadfield matrix as input to the exact low-resolution brain electromagnetic tomography algorithm (Pascual-Marqui et al 2011), in order to reconstruct neural signals in the grey matter. The signals were then transformed to Montreal Neurological Institute (MNI) space by using a non-rigid deformation, calculated with SPM12 (Ashburner et al 2014) from the T1-weighted individual image.

For each voxel in MNI space, we then performed a time-frequency decomposition of the estimated time-course by using short-time Fourier transform with 1 s sliding window and 0.1 s steps, and based on that we performed an ERD/ERS analysis. To this end, the spectrogram was epoched using a time window between −1 and +4 s with respect to the movement onset, and averaged. The onset of each movement block was defined as the maximum of the first derivative of the absolute acceleration, averaged across all the corresponding channels. Finally, the ERD/ERS (Pfurtscheller and Da Silva 1999) was quantified using the following formula:

$$\begin{equation}{\text{ERD/ERS}}\left( {f,t} \right) = \frac{{P\left( {{\,f},t} \right) - {P_{\text{B}}}\left( {f\,} \right)}}{{{P_{\text{B}}}\left( {f\,} \right)}} \times 100\% \end{equation} \tag{ 4 }$$

where $P\left( {{\,f},t} \right)$ is the power at any given frequency and time, and ${P_{\text{B}}}\left( {f\,} \right)$ is the average power in a baseline period, defined between −1 and −0.1 s with respect to the movement onset (Zhao et al 2019). ERD/ERS values, which were initially calculated for each voxel, were then extracted from specific regions of interests (ROIs): bilateral primary motor cortex (M1, MNI coordinates: [±44, −17, 49]) and thalamus (MNI coordinates: [±9, −12, 10]). These ROIs were selected as they are relevant for the motor task under investigation (Fukuyama et al 1997, Zhao et al 2019). The thalamus was also included, as it has been shown that hdEEG may permit the reliable reconstruction of neural activity in both cortical and subcortical regions (Zhao et al 2019, Samogin et al 2020). Furthermore, ERD/ERS values were quantified by averaging, either for the ROIs or for individual voxels, across latencies between +0.1 and +4 s, and across frequencies in the alpha [8–13 Hz] and beta [13–30 Hz] bands (Zhao et al 2019, Guarnieri et al 2021). ERD/ERS values obtained across individual voxels were used to produce ERD/ERS maps.

For each artifact removal approach, we calculated the spatial correlation between the ERD/ERS maps in the walking conditions, and a template ERD/ERS map of foot movement, separately for the alpha and beta bands. This template map was bilateral, and obtained by averaging the ERD/ERS map of the right foot movement in a sitting condition (Zhao et al 2019) with the same map flipped along the left-right axis. This template map was used as reference, as it was obtained under controlled movement conditions, and was considered to be largely unaffected by artifacts.

2.2.4. Statistical analyses

We used an analysis of variance (ANOVA) to examine the performance of different processing method (multi-step BSS, FastICA-defl, Infomax, no processing) on the noise and signal levels (from auditory hdEEG data), ERD intensities and spatial correlations (from motor hdEEG data). The correlations values were transformed to z-values using the Fisher transformation, to improve the data normality. A two-way ANOVA was used to test the impact of the computational approach and the experimental condition on ERP noise and signal level, respectively. A three-way ANOVA was used to test the impact of the computational approach, the frequency band and the experimental condition on ERD/ERS intensities and spatial correlations, respectively. For post-hoc comparisons between computational approaches, we used two-side Wilcoxon signed rank tests. The significance level was corrected for multiple comparison using the false discovery rate (FDR) method. Group-level analysis were conducted by averaging values across participants.

3.1. Definition of multi-step BSS settings

Our initial analyses were directed toward the choice of a specific BSS method for the three processing steps. Based on the maximum absolute correlation with reference signals, we selected FastICA-defl for the first step (ocular artifact attenuation), FastICA-symm for the second step (movement artifact attenuation), and IVA-G for the third step (myogenic artifact attenuation) (figure 2). We then defined the thresholds that permitted to optimally separate artifactual and non-artifactual components in each step. The window size for kurtosis calculation (first step) was set to 5 s and the one for sample entropy (second step) to 20 s. Furthermore, the cut-off frequency for the definition of the high-frequency band (third step) was set to 30 Hz. The classification thresholds defined for the three processing steps were equal to 12 (maximum kurtosis), 0.78 (mean sample entropy) and 0.53 (relatively high-frequency power), respectively. Based on the defined BSS methods and classification parameters, we estimated to have 92.07%, 87.28% and 93.54% classification accuracy for the processing steps, respectively (figure 3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2. Artifact attenuation performance

When using the proposed BSS approach, we expected to observe an effective attenuation of ocular, movement and myogenic artifacts. For qualitative analysis, we displayed EEG traces for each processing step (figure 4), focusing on the data segments in which the artifacts were prominent before correction. We observed that the intensity of the artifacts after processing was reduced to the point that no clear residual could be detected. To corroborate this result, we then conducted a quantitative analysis based on the correlation between EEG signals and reference artifactual signals. For all the task conditions, the maximum correlation between EEG and EOG signals showed a significant decrease ($p < 0.001$) following ocular artifact correction (figure 4(a)). Also, the maximum correlation between EEG and limb velocity signals showed a significant decrease ($p = 0.002$ and $p < 0.001$ for slow and fast ground walking conditions respectively) following movement artifact correction (figure 4(b)). Finally, the maximum correlation between EEG and EMG signals showed a significant decrease ($p < 0.001$ for both slow and fast ground walking conditions respectively) following myogenic artifact correction (figure 4(c)).

Zoom In Zoom Out Reset image size

3.3. Analysis of ERP results

Next, we compared the multi-step BSS approach for EEG artifact attenuation to other solutions by using auditory-evoked ERP data. Visual examination of ERPs confirmed the presence of a strong response to auditory stimulation at the latency of about 100 ms, also known as N1 response. This response was generally most prominent over medial central regions. Different levels of noise in the pre-stimulus interval were present depending on the specific data processing used, although they were relatively difficult to appreciate by the naked eye (figure 5).

Zoom In Zoom Out Reset image size

Quantitative analyses evidenced differences among artifact attenuation methods (figure 6). A two-way ANOVA indicated significant differences in noise level across methods ($p < 0.001$) and conditions ($p = 0.023$), with no significant interaction ($p = 0.916$) between method and condition (table S1). In the stance condition, the noise level with multi-step BSS was lower compared to the unprocessed data (${p_{{\text{FDR}}}} = 0.002$) and other approaches (${p_{{\text{FDR}}}} = 0.005$ for FastICA-defl, and ${p_{{\text{FDR}}}} = 0.046{ }$ for Infomax, figure 6(a)). Even stronger differences (${p_{{\text{FDR}}}} < 0.001{ }$ for multi-step BSS with respect to other methods) were found also for slow and fast ground walking conditions, respectively. No significant differences were found in signal level with different approaches ($p = 0.178$) and conditions ($p = 0.075$), as revealed by a two-way ANOVA (table S2).

Zoom In Zoom Out Reset image size

3.4. Analysis of ERD/ERS results

Further evidence for the effectiveness of multi-step BSS was gathered by investigating ERD/ERS during slow- and fast-walking conditions. We observed indeed an ERD in alpha and beta bands, both for M1 and thalamus, by using multi-step BSS. Before artifact removal, instead, no evident modulation of band-limited neural power was observed in the fast-walking condition, and only weak effects could be seen in the slow-walking condition. The ERD obtained for M1 and thalamus with FastICA-defl and Infomax, respectively, were detectable, but less strong than that discerned with multi-step BSS (figure 7).

Zoom In Zoom Out Reset image size

The intensities of the alpha and beta ERD from individual participants and ROIs were extracted and analyzed at the group level using a three-way ANOVA (figure 8 and table S3). This analysis revealed significant differences in ERD across artifact attenuation methods ($p = 0.012$) and between movement conditions ($p < 0.001$). No significant differences were found between alpha and beta bands ($p = 0.072$). In the slow-walking condition, paired comparisons showed significantly stronger alpha and beta ERDs with multi-step BSS than without artifact correction (${p_{{\text{FDR}}}} = 0.010$ for both bands) and with FastICA-defl (${p_{{\text{FDR}}}} = 0.010$ and ${p_{{\text{FDR}}}} = 0.043$, respectively). In the fast-walking condition, the data processed with multi-step BSS had stronger alpha and beta ERD than the unprocessed data (${p_{{\text{FDR}}}} = 0.003$ and ${p_{{\text{FDR}}}} < 0.001$, respectively), as well as than the data cleaned using FastICA-defl (${p_{{\text{FDR}}}} = 0.004$ and ${p_{{\text{FDR}}}} < 0.001$, respectively) and Infomax (${p_{{\text{FDR}}}} = 0.007$ and ${p_{{\text{FDR}}}} < 0.001$, respectively).

Zoom In Zoom Out Reset image size

As a final analysis, we examined the ERD/ERS source maps produced using different artifact removal methods. The alpha ERD map was relatively widespread, with strongest values over central-parietal regions. Overall, ERD intensity was larger with multi-step BSS than FastICA-defl and Infomax, both during slow- and fast-walking (figure 9). When examining the beta ERD map, we found that multi-step BSS allowed identifying the premotor area, the supplementary motor area, and the dorsal part of the primary motor cortex (figure 9(b)). For alpha and beta bands separately, we quantified the spatial correlation of the ERD maps from each individual, against a motor-related template map (figure 10). A three-way ANOVA (table S4) showed significant differences in spatial correlation values for method ($p < 0.001$) and band ($p < 0.001$). No significant difference was found between task conditions ($p = 0.798$). For all conditions and bands, multi-step BSS yielded higher spatial correlations. In the slow-walking condition, we did not find significant differences as compared to other approaches. In the fast-walking condition, higher correlations were observed for both bands in comparison to Infomax (${p_{{\text{FDR}}}} = 0.010$ for the alpha and ${p_{{\text{FDR}}}} = 0.015$ for the beta band). No significant difference was observed as compared to FastICA-defl. Similar results were found when examining the spatial similarity between the group-level ERD maps and the motor-related template maps in alpha and beta bands, respectively (figure S3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The study of neural dynamics in naturalistic conditions, such as for instance during walking, is gaining an increasing interest by the neuroscientific community. To support this line of research, novel mobile brain/body imaging applications based on EEG are being developed (Gramann et al 2011, Gwin et al 2011, Seeber et al 2014, 2015, Wagner et al 2016, Artoni et al 2017, Bradford et al 2019, Nakagome et al 2020). However, the interpretation of the results emerging from mobile hdEEG recordings is strongly challenged by the presence of movement and myogenic artifacts (Nathan and Contreras-Vidal 2016, Richer et al 2020). The main goal of this study was to develop a multi-step BSS solution for attenuating ocular, movement and myogenic artifacts, which are present in mobile hdEEG data (Fitzgibbon et al 2007, Nathan and Contreras-Vidal 2016). The order with which those artifacts were attenuated was defined on the basis of their relative power (Urigüen and Garcia-Zapirain 2015, Richer et al 2020). For each step, we defined the BSS method to be used, among PCA, CCA, ICA (FastICA, Infomax and JADE) and IVA. The selected method was the one with the best decomposition performance, measured in terms of maximum absolute correlation between reference artifactual signals and component signals. The results we obtained in terms of BSS method selection are consistent with previous studies suggesting that ICA can be effectively used for ocular (Guarnieri et al 2018, Barban et al 2021) and movement artifacts (Snyder et al 2015), as well as IVA outperforms other solutions in the case of myogenic artifacts (Chen et al 2017, Barban et al 2021). After the selection of the BSS method, we defined proper thresholds for the classification of the components, using specific parameters depending on the kind of artifact to be removed. In particular, we measured kurtosis and sample entropy across consecutive data segments to extract time-dependent features of ocular and movement artifacts. As for the classification of the myogenic artifact, we quantified the ratio between the power in the gamma band and the total band, considering that the power of neural signals decreases as the frequency increases (Daly et al 2012). It should also be noted that, since the artifact classification for each step is based on a single parameter, our multi-step approach can be easily adapted, if needed, by modifying the classification thresholds. In sum, we observed that the multi-step BSS approach is superior to single-step BSS approaches, not only because it allows optimal separation of different kinds of artifacts from neural activity, but also because each artifactual source can be easily detected based on an optimized parameter-based thresholding.

The effectiveness of artifact attenuation by multi-step BSS was demonstrated in first instance by the low values of correlation between artifact-corrected hdEEG signals and reference artifactual signals. Specifically, the median value of maximum absolute correlation between EEG and EOG signals dropped—on average across conditions—from 0.92 to 0.06 after artifact attenuation. Furthermore, the maximum absolute correlation between EEG and limb velocity signals was reduced from 0.26 to 0.08. A similar reduction after artifact attenuation was also observed for the maximum absolute correlation between the envelopes of EEG and neck EMG signals, which went from 0.29 to 0.01.

The performance comparison of multi-step BSS with respect to single-step BSS approaches (FastICA-defl and Infomax) using ERP and ERD/ERS analyses corroborated the results described above. Both FastICA-defl and Infomax were chosen as benchmarks, as they are widely used for artifact removal, and are paired with effective solutions of automated classification of artifactual sources (Liu et al 2017, Pion-Tonachini et al 2019). We did not directly test the performance of two-step BSS procedures in this study. However, our analyses suggested that ocular, movement and myogenic artifacts can be removed more effectively with different BSS methods, and with different artifact classification features. We can then infer that using a single BSS method and a single classification feature for two of the three artifacts considered in this study, as for instance movement and myogenic artifacts, may lead to suboptimal results.

The use of automated classification tools not only facilitated our computational analyses, but also allow us to obtain reproducible results, unaffected by manual classification (Mantini et al 2008, Pion-Tonachini et al 2019). Moreover, by using auditory ERP data, we quantified residual noise levels present in the pre-stimulus period, as we previously did in a number of studies (Guarnieri et al 2018). For all the auditory conditions, multi-step BSS showed the lowest values, with significant differences as compared to other methods. Examining the ERP data, we also noticed that Infomax has slightly stronger N1 peak than with other approaches. This method seemed to produce different topography of N1 responses compared to previous reports (Giard et al 1988, Liebenthal et al 2003, Zink et al 2016). We speculate that the artifact attenuation method based on Infomax may introduce low-frequency components in the reconstructed EEG signal, thereby altering the N1 intensity over the scalp. Such distortions were not observed with FastICA-defl and—most importantly—with multi-step BSS.

By examining ERD/ERS in alternated walking-stance conditions, we evaluated the modulations of neural power at the level of specific ROIs as well as of the whole brain. This was done using advanced techniques for neural activity reconstruction from hdEEG data, based on a realistic individual head model (Liu et al 2017, Taberna et al 2021a, 2021b). Overall, the ERD/ERS analysis evidenced that the attenuation of artifacts with multi-step BSS permit to recover movement-related neural dynamics in naturalistic conditions. Differences between multi-step BSS and other solutions under evaluation, namely FastICA-defl and Infomax, were especially observed in the fast-walking condition, which is characterized by the presence of stronger artifacts in hdEEG data. Importantly, ERD/ERS maps obtained from hdEEG data processed using multi-step BSS showed, particularly for the beta band, the involvement of a widely distributed cortical network of gait-related regions (Seeber et al 2014, Knaepen et al 2015, Hinton et al 2019). The use of other artifact attenuation methods, in turn, yielded ERS/ERD maps, whose primary regions cannot be easily related to motor control processes (Zhao et al 2019). Overall, these results suggest that multi-step BSS may be essential to reliably assess neural dynamics in mobile hdEEG experiments. Interestingly, by using mobile hdEEG we could observe differences in ERD intensity between slow- and fast-walking conditions. These differences might be directly linked to neural processes supporting gait control (Bulea et al 2015); alternatively, they might be due to a different level of artifactual residuals in mobile hdEEG signals for slow- and fast-walking conditions, respectively.

A number of limitations of our study should be mentioned. Firstly, we did not test the generalizability of our approach for a wide range of tasks, including more complex motor tasks. Future studies are therefore warranted to ensure that the selected BSS methods and parameters for each step may yield satisfactory artifact attenuation for different EEG datasets. Secondly, we did not assess any approach for the attenuation of the EEG artifacts that are present during functional magnetic resonance imaging, such as ballistocardiographic artifacts (Marino et al 2018). Our approach can in principle be extended to such applications in the future. Lastly, we used a single parameter for the classification of artifactual sources following each BSS step. Notably, we posit that the performance of component classification, hence the general effectiveness of our approach, may be further improved by applying machine learning tools, as for example deep learning (Zhang et al 2021a, 2021b).

We have introduced a novel multi-step BSS approach for attenuating ocular, movement and myogenic artifacts, which are typically present in mobile hdEEG data. For each step, one of these three artifact categories is specifically addressed, using a specific BSS method and an optimized artifact classification parameter. Our results showed that the use of the novel approach yielded lower noise in the hdEEG data, and permitted to retrieve stronger task-related modulations of neural activity than alternative solutions. Our computational solution supports a wider use of hdEEG-based source imaging in movement and rehabilitation studies, and may contribute to the further development of mobile brain/body imaging applications (Jungnickel et al 2019).

The work was supported by the KU Leuven Special Research Fund (Grant C16/15/070), the Research Foundation Flanders (FWO) (Grants G0F76.16N, G0936.16N and EOS.30446199 to D M, and fellowship 1104520N to R G), the Chinese Scholarship Council (Scholarship 201708620182 to M Z), and the Italian Ministry of Health (Grants RF-2018-12366899 and SG-2018-12368232).

The data that support the findings of this study are available upon reasonable request from the authors.

The authors declare no competing interests.