For this study, we relied on a cybersickness dataset described in Li et al. (2021), which followed an experimental protocol specifically designed to induce cybersickness. Data was collected from 20 participants who were immersed in two VR scenarios: a tunnel travel simulation and a roller coaster ride. The ‘Tunnel’ session emphasized strong linear vection with minimal rotational motion to induce mild to moderate cybersickness. Participants who did not reach a cybersickness threshold of 11/20 on the Fast Motion Sickness (FMS) scale after 10 repetitions (approximately 10 min) proceeded to the ‘Roller Coaster’ session. This session included similar linear vection but added stronger rotational motion to intensify symptoms, ensuring that all participants reached sufficient cybersickness levels to enable subsequent analyses. The rationale behind these tasks was to account for individual differences in susceptibility while ensuring sensory conflict through visual motion without movements.

Each session lasted up to 10 min, with the duration tailored to individual participant responses, as indicated by their cybersickness levels. Cybersickness severity was systematically measured using the 20-point FMS questionnaire, administered at the end of each minute, for as long as the participants could manage their symptoms until reaching an FMS score of 11. Some participants did not experience cybersickness and reported FMS values not exceeding 4 (out of 20). Specifically, subjects 1, 3, 10 and 15 in the Tunnel session and subjects 1, 3, 10 and 16 in the Roller coaster session. In addition, data from 10 subjects in the Roller Coaster session was not available since these subjects did not continue to the second part of the experiment due to severe discomfort after the first video. In the roller coaster session only the data from subjects 7, 8, 11, 13, 15 and 20 were included in the analysis.

The study recorded multi-channel EEG, 1-channel Electrooculography (EOG), and head movements via a 3-axis accelerometer. The seven EEG channel locations used in this study included Fz, Pz, P3, P4, Cz, CP5, and CP6, recording with a sampling rate of 500 Hz. The EOG channel, located on the left lower eyelid, was used to monitor eye movements and eye blinks. Both the ground and reference electrodes were connected together and attached to the right earlobe using an ear clip. The interested reader is referred to Li et al. (2021) for complete details about the data collection process and the dataset used.

Pre-processing is crucial in EEG research to remove unwanted physiological (e.g., heart beats, muscle/eye/head movements) and environmental (e.g., powerline, temperature, humidity) interferences (Tanner et al., 2016; Bigdely-Shamlo et al., 2020). With such pre-processing algorithms, a trade-off must always be achieved between the aggressiveness of the artifact removal step and the potential (unwanted) removal of useful discriminatory information from the collected signals (Delorme, 2023; Robbins et al., 2020). For example, AAR algorithms based on Independent Component Analysis (ICA) are susceptible to causing phase distortion (Thatcher et al., 2020) and/or loss of useful information (Sadiya et al., 2021; Klug and Gramann, 2021). Similarly, rejecting channels and segments of the signal contaminated by artifacts introduces information loss and discontinuities (Sadiya et al., 2021; Huberty et al., 2024).

With cybersickness prediction, however, artifact removal may play a crucial role, as the information that is often removed, such as head movements, eye movements, and heart beats, may indeed help in the prediction process. As such, here we explore three different levels of AAR, from light to aggressive, to gauge the effect this has on cybersickness characterization, as well as on the insights generated by the top features later used for the prediction task. This, in fact, has been recently emphasized in Delorme (2023) where a reduction in the amount of pre-processing was suggested to prevent the loss of relevant neural information of event related potentials. In addition, authors in Thatcher et al. (2020) demonstrate that AAR based on ICA result in phase distortion that could compromise connectivity analysis. The three different pre-processing techniques explored herein are detailed next.

Figure 2 illustrates the general testing methodology employed in this experiment. With the dataset used, cybersickness FMS ratings were provided every minute.

FMS scores were recorded by participants at the end of each 1-min interval, while the EEG data used in regression analysis were taken from the same 1-min segment. This means that EEG trials could correspond to any moment within the interval-potentially near the beginning or end of the minute-relative to the reported FMS score. This approach may introduce a small temporal mismatch between the EEG features and the moment the FMS score reflects. As such, in this study we are measuring the cybersickness level (given by the FMS score) on a per-minute basis. Future studies could take this further and use the EEG trials to predict cybersickness levels many minutes ahead.

To augment the amount of data available for training of the classifiers, a bootstrapping method was applied where a random selection of epochs within a certain 1-min EEG segment under the same FMS rating were used. The statistical functionals, such as mean, standard deviation, and skewness were then computed from the resulting subset. This process can be repeated multiple times where after each iteration a distinct set is created. The number of epochs taken per 1-min segment and the number of times this bootstrap sampling was done was empirically set to 75% of total epoch number within the minute and 10 repetitions, respectively.

The data augmentation process produced a total of 712 samples from the collected data, as compared to the 7,035 features extracted for analysis. Feature selection methods were implemented to mitigate the high sample-to-feature ratio, ensuring a more robust and interpretable regression model. The tuning stage involved the use of a validation set to estimate the best subset of features, as depicted in Figures 2A, B.

To avoid data leakage between training and test samples, a Stratified-group K-fold methodology was used to restrict groups of samples from the same minute segment to a unique data subset (i.e., train or test). This step corresponds to item (d) in Figure 2. Train and test set were generated from the development set which correspond to a portion of 67% of samples reserved exclusively for training and testing. Hyperparameters of the support vector regression model were optimized using the validation set, as represented by Figure 2C, during the tuning stage.

After extracting features from each EEG recording session, outlier detection was conducted independently on each feature using a threshold of three times the median absolute deviation. Detected outliers, along with near infinite, near zero, and missing samples (Not-A-Number), were corrected via linear interpolation of time contiguous samples. Lastly, FMS ratings were normalized to a scale ranging from 0 to 1 to facilitate the use of the root-mean-square error (RMSE) as a figure-of-merit to gauge system accuracy.

Given that the used dataset has a low sample-to-feature ratio, feature selection is needed to reduce the number of features to a manageable number to prevent overfitting. Additionally, to mitigate the risk of data leakage, features selection is performed on a validation subset (set as 33% samples of the dataset) of the data using the “StratifiedGroupKFold” module within scikit-learn (Pedregosa et al., 2011) prior to the train/test split procedure. This approach ensures that the validation samples come from different time instances than those reserved for testing, keeping the test set completely unseen.

Here, a simple feature ranking method was utilized based on Spearman correlation. For each feature, the correlation was calculated between samples and their corresponding FMS rating. This process is represented by step (a) in Figure 2. The correlation score attributed to each feature is used in ablation study II presented in Section 3.2.

The rank ordering was based firstly by statistical significance (Holm-Bonferroni corrected with significance set at a threshold of 0.05), followed by absolute Spearman correlation coefficient values. After this ranking, the top 100 features were selected (empirically) for further analysis. To gauge the impact of including additional features on cybersickness measurement, top-features were added one by one and regression performance on the validation set was assessed via stratified 5-fold cross-validation. This step corresponds to item (b) of the Tuning stage presented in Figure 2; this step further described in the ablation study III in Section 3.3. The final number of features to be used is the one that leads to the smallest RMSE value obtained on the validation set after Gaussian smoothing.

A Support Vector Regressor (SVR) was selected as the regression algorithm, motivated by its robustness to overfitting, a characteristic particularly advantageous for high dimensional datasets. Hyper-parameter optimization for the SVR model was carried out through an exhaustive 3-fold cross-validation grid search, exploring combinations of 3 × 4 × 20x20 for kernel types (‘poly’, ‘rbf’, and ‘sigmoid’), polynomial degrees (2 through 5), and the parameters C and gamma. The search covered logarithmic scales from 1 0 − 4 to 1 for C, and from 1 0 − 4 to 1 0 2 for gamma, with ‘auto’ and ‘scale’ additionally evaluated as specific gamma parameters. This step is represented by block (c) in Figure 2. Prior to training the SVR, data was normalized subtracting the mean and scaling to unit variance to address SVR sensitivity to scaling. Bootstrapping was done 100 times on the 5-fold cross-validation task to allow for significance calculation of the obtained RMSE values relative to chance.

To better understand the impact of head movements and eye movements on certain EEG frequency bands, Figures 3, 4 depict (top plots) the power spectral representation of the accelerometer signals (averaged across the x, y, and z-axes) and the Signal-to-Noise Ratio (SNR) representation (bottom plot), for each of these two signal modalities, respectively. SNR is computed for each frequency bin as the ratio between the power of the bin’s central frequency and the mean power of neighboring frequencies excluding immediately adjacent ones (Meigen and Bach, 1999).

As can be seen from Figure 3, notable power peaks occur at 1.7 Hz, 3.3 Hz, and 4.9 Hz, likely harmonics of the head movements. These findings suggest that head movement information can likely be a confound for EEG features based on delta and theta bands. Similarly, from Figure 4 ocular activity exhibits a 1/f spectral roll-off, suggesting that eye movements may also affect lower frequency bands, specifically the delta and theta bands.

Next, we wish to explore the potential of the three different AAR algorithms at removing these two artifact sources. To this end, we utilize the mutual information computed from the processed EEG signal and the temporal series of the EOG and the averaged tri-axes accelerometer signals. A higher mutual information will signal that the artifacts have not been removed and the EEG signal still contains substantial artifactual information embedded in the signal. In turn, lower mutual information values suggest that the artifacts were correctly removed. Figure 5 presents the computed mutual information (averaged over the head and eye movement signals) for each of the seven EEG channels. As can be seen, minimal processing showed the highest mutual information, as expected. The RegICA method removed some of the artifacts, while the more complex ASR + ICLabel combination resulted in the least amount of mutual information, suggesting the most aggressive removal of artifacts.

Lastly, Figure 6 depicts the EEG Power Spectral Density averaged across channels processed by the Minimal and the ASR + ICLabel pipelines. As can be seen, the 1.7 Hz, 2.5, and 3.2 Hz peaks, likely due to head movements, can be seen with the former, but not the latter. The spectral density above 6 Hz, in turn, closely matches in both scenarios.

Comparison of the spectra of the two pipelines, Minimal and ASR + ICLabel, suggests that the components identified as ocular and motion artifacts, in Figures 3, 4 respectively, were eliminated by the ASR + ICLabel pipeline, but not by the Minimal pipeline. This observation in the spectral domain is consistent with the results obtained from the mutual information measurements.

Despite the effective removal of artifacts, AARs can distort the signals. Therefore, we investigated their impact on the useful information related to cybersickness. We compared different combinations of epoch sizes and each AARs for the predictive power of each type of feature. As a figure-of-merit, we employed the correlation measures described in Section 2.6, which were calculated between the samples and the associated FMS rating independently for each feature.

Figure 7 presents the Spearman correlation between different features and the FMS scores in the training set. The correlation map shows results for the different AAR methods, feature types, and epoch sizes. As can be seen, PSD and AMP features showed the highest correlations with the Minimal processing pipeline with the 8 and 20 s epochs for AMP features and 8-s epochs for PSDs. Moreover, the band and channel coupling measures showed mild to moderate correlation with FMS across various pre-processing scenarios, suggesting their potential complementarity for cybersickness characterization.

Overall, the RegICA and ASR + ICLabel AAR methods showed lower correlation scores with FMS ratings across most epoch sizes and feature types, whereas the Minimal pre-processing showed the highest. These findings suggest that, indeed, retaining artifacts related to blinks, eye movements, and head movements can be useful to characterize the multimodal effects of cybersickness.

Next, we explore the selection of the optimal number of features. Figure 8 shows the coefficient of determination ( R 2 ) scores as a function of number of features for the case of ASR + ICLabel pre-processing and 5-s epochs. As can be seen, the R 2 score improves with the inclusion of more features until a plateau is reached, which starts with a R 2 score of 0.467 at around 71 selected features. The optimal feature number is selected for each pre-processing and epoch size combination; the exact number of features selected for each of these combinations is shown in Table 1. These top features are discussed next.

To gain greater insight into cybersickness characterization, two analyses were performed based on the interpretation of the most relevant features. The first analysis explores the neural patterns indicative of cybersickness, while the second investigates the multi-modal potential of cybersickness detection with minimal EEG pre-processing.

We employed a strategy to summarize the multiple lists of top features generated for each combination of the AAR algorithm and epoch size. For each AAR, a unified list was created by gathering all features that appeared in at least one of the lists across various epoch sizes. A common list is then defined as the intersection of two unified lists, meaning that a feature must be present in both lists to be included. Therefore, the analysis in the discussion section is based on common features between the two AARs, regardless of the epoch size.

Table 3 lists the final RMSE achieved between the real (normalized) FMS ratings and the predicted ones for each pre-processing strategy and epoch size. These results reflect the model’s performance on unseen data, as the training and testing folds used in cross-validation are strictly non-overlapping. To gauge the significance of the results in the table over chance, a random classifier is also trained where the FMS ratings, which serve as regression labels, are randomized. As can be seen, across all tested configurations, performance significantly exceeded this chance level. Table 4 summarizes the main relevant features used for cybersickness level measurement.

Overall, the ‘Minimal’ pre-processing approach achieved the lowest RMSE values with top features computed from 8-s epochs, achieving an RMSE of 0.234, significantly lower than chance-level (0.327), corroborating the importance of the multi-modal aspect achieved by leaving artifacts untouched. In contrast, the RegICA and ASR + ICLabel pre-processing approaches achieved their lowest RMSE values at higher epoch sizes. As these methods tend to remove artifacts, the regressors are relying solely on neural signatures for cybersickness prediction. In such cases, the longer duration epochs are needed to capture the neural underpinnings of cybersickness, which may be present in lower-frequency activity.

In EEG studies, the presence of artifacts can negatively influence the validity and interpretation of results (Zhou et al., 2023; Lopes et al., 2023; Bennett et al., 2021), as artifacts can take on various forms (Delorme, 2023; Michel and Brunet, 2019; Jin et al., 2023) and overlap with neurological signals in spatial, temporal and spectral domains (Thompson et al., 2008; Chiarion et al., 2023). While this makes the use of pre-processing an essential step in EEG signal analysis, there is no universally accepted method. Algorithms are developed to target certain types of artifacts [e.g. (Robbins et al., 2020; Yu et al., 2022)] and each pre-processing method comes with its own set of advantages and limitations. As such, researchers must trade-off between effective artifact removal and minimal signal distortion (Bouazizi and Ltifi, 2024; Tajmirriahi et al., 2022; Bao et al., 2022).

Recently, the extensive use of AAR has been contested due to the excessive data distortions and destruction they might cause to the data (Delorme, 2023; Chiarion et al., 2023). For instance, ICA has been shown to produce distortions due to the non-stationary nature of EEG signals (Chiarion et al., 2023; Bouazizi and Ltifi, 2024). Additionally, it may mix neurological and non-neurological sources within the same components when the number of channels is much smaller than the number of sources (Thompson et al., 2008; Sun and Mou, 2023; Chaddad et al., 2023), which can often be the case in VR applications.

As with most EEG studies, our results also indicate that non-neurological activity is recorded along with neurological activity. More specifically, we show that ocular activity and head movements are observed in the lower frequencies of the EEG spectrum, such as delta and theta, and these were selected as top-features in the Minimal processing case. In turn, other features from these regions, but measuring temporal dynamics information, could be indicative of regulatory and coordinating roles in cognitive processes.

Unlike other studies, however, for cybersickness prediction, the presence of physiological artifacts can serve as useful sources of multimodal information to predict cybersickness levels (Islam et al., 2021; Dennison et al., 2016; Islam et al., 2022; Chang et al., 2021; Shimada et al., 2023b; Jeong et al., 2022). The obtained findings showed that Minimal processing, indeed, resulted in the lowest cybersickness level prediction error. Therefore, if one is interested in maximizing cybersickness prediction accuracy, performing minimal pre-processing may be the way to go. It is suggested, however, that careful scrutiny of the top-used features be evaluated, to ensure that indeed the developed models are also relying on relevant neural data, and not just physiological artifacts. Furthermore, the main sources of artifacts should be precisely identified when implementing an EEG-embedded multimodal system, as not all artifacts are useful for predictions and may negatively impact outcomes.

This research presents some limitations. First, the use of a non-linear SVR kernel made it difficult to analyze feature ranking by inherent regularization, which is possible only with linear kernels. As a result, we relied on Spearman correlation for feature ranking, which does not consider the interaction between features. Future work may explore the use of more advanced feature selection methods that take such interactions into account, such as the mRMR (Minimum Redundancy Maximum Relevance) (Ding and Peng, 2005).

The findings of this study apply primarily to open-loop VR systems, where users have limited control over motion stimuli. EEG patterns observed under these conditions may not directly translate to closed-loop systems, where active control relies on anticipatory mechanisms, as shown in research comparing drivers and passengers (Dong et al., 2011; Rolnick and Lubow, 1991), in contrast to the compensatory responses observed in open-loop scenarios. While the study provides valuable insights into the neural correlates of cybersickness, its findings are primarily applicable to open-loop VR systems. Given that closed-loop systems are more common among general VR users, further research is needed to validate the relevance of these markers in user-controlled environments.

Moreover, here we relied on the efficiency metrics for pre-processing algorithms based on the mutual information difference with head accelerometer and one EOG time serie. Other artifacts may contaminate the EEG signals, such as muscle movements, poor electrode-skin contact, and sweat, to name few. Continuous measurement of such signals are more challenging, but alternate modalities, such as galvanic skin response, facial electromyography (EMG), or impedance measurements may provide additional cues for future works. Emerging VR headsets are already including such modalities directly into the head-mounted display, thus facilitating this type of analyses in the future. In this study, we relied on a publicly available dataset, which primarily included artifacts from ocular and head movements. While these are common in VR scenarios, future work should validate the impact of minimal preprocessing under a broader range of artifact types to ensure the generalizability of these findings.

Lastly, the proposed work relied on publicly available data, which was collected from a limited number of participants, and not all of them completed the full task. This limitation led to variability across subjects, particularly at higher FMS rating levels. Additionally, the dataset’s gender composition (31 females, 9 males) introduces a potential bias, as the findings may predominantly reflect female responses, consistent with studies reporting higher susceptibility to motion sickness among females (Koslucher et al., 2015; Munafo et al., 2017). Furthermore, while this study focused on the sensory conflict theory, other theories, such as the ecological theory of motion sickness and postural instability (Stanney et al., 2020; Riccio and Stoffregen, 1991), emphasize individual differences and environmental interactions that were beyond the scope of this work. Future research could investigate these complementary perspectives and employ advanced cross-subject generalization tools (e.g. (Albuquerque et al., 2022)) to address these limitations.