Subjective Fear in Virtual Reality: A Linear Mixed-Effects Analysis of Skin Conductance

The investigation of the physiological and pathological processes involved in fear perception is complicated due to the difficulties in reliably eliciting and measuring the complex construct of fear. This study proposes a novel approach to induce and measure subjective fear and its physiological correlates combining virtual reality (VR) with a mixed-effects model based on skin conductance (SC). Specifically, we developed a new VR scenario applying specific guidelines derived from horror movies and video games. Such a VR environment was used to induce fear in eighteen volunteers in an experimental protocol, including two relaxation scenarios and a neutral virtual environment. The SC signal was acquired throughout the experiment, and after each virtual scenario, the emotional state and fear perception level were assessed using psychometric scales. We statistically evaluated the greatest sympathetic activation induced by the fearful scenario compared to the others, showing significant results for most SC-derived features. Finally, we developed a rigorous mixed-effects model to explain the perceived fear as a function of the SC features. Model-fitting results showed a significant relationship between the fear perception scores and a combination of features extracted from both fast- and slow-varying SC components, proposing a novel solution for a more objective fear assessment.


INTRODUCTION
F EAR is an emotion reported by persons who feel afraid for the actual or potential presence of danger, pain, or harm [1], [2]. Following the exposure to fearful stimuli, psychophysiological response patterns can be recorded. This led to infer a conscious feeling of fear based on the detection of its mere physiological correlates -even when the subjective emotion has not been experienced (e.g., in response to subliminal phobic stimuli; see [3] for a comprehensive review). To overcome this ambiguity, the subjective emotional state consciously reportable (fear) has been proposed to be distinguished from the typically -but not inevitablyassociated physiological correlates [2]. The subjective feeling of fear can be directly assessed only through the subject introspection [2], however self-report measures are susceptible to potential issues (see [4] for a critical review of one of the most widespread self-report tools). For example, people may lie about their feeling for a variety of reasons (e.g., shame), or they could just lack the introspective abilities needed to recognize their feelings, since these are regulated mainly (or, at least, also) out of conscious control [5]. Such unreliability led researchers to devise more objective and replicable measures of the emotional state by analyzing correlates of sympathetic nervous system (SNS) activity, the increase of which induces the bodily changes necessary to cope with environmental threats [6], [7]. One of the most studied physiological correlates of SNS activity is the skin conductance (SC) of the hands, which has been considered an ideal way to estimate the human arousal level [8]. Indeed, due to its direct and exclusive control performed by the SNS, the SC provides a relatively direct and undiluted representation of sympathetic activity.
Despite the wide use of SC measurements in the scientific literature studying emotional reactions [9], an extensive and rigorous investigation of the SC components mainly correlated to subjective fear is still missing. Nevertheless, including such information in a computational model able to continuously and objectively estimate the perceived fear could have important implications for phobia clinical therapies. First, a more objective measure of fear could allow a better assessment of symptom severity and complexitycurrently based mostly on self-reports susceptible to psychological variables. Second, such a model would allow a real-time assessment and adjustment of the phobic stimulus in exposure therapies, improving both the efficacy and acceptability of the treatment, thanks to a more gradual and personalized exposure.
One of the biggest difficulties in this research field concerns the possibility of eliciting specific emotions in an effective and realistic -but also replicable and controllable -way, and in analyzing complex variance structures often associated with such kind of experimental data. In fact, the most widespread scientific tools to elicit emotional reactions are based on a non-immersive, passive media exposure. Standardized sets such as the International Affective Picture System (IAPS), the International Affective Digitized Sounds (IADS), or the film-clip paradigm have often been unable to induce discrete emotions such as fear, anger, and disgust, that are not easily separable from each other [10], [11]. In addition, previous studies have claimed that passive exposure to images or sounds can hardly reach a level of engagement sufficient to selectively elicit emotions in a controlled and reproducible way [12]: seeing a picture of someone holding a knife is very different from being personally menaced by an armed person. This gap between setting's reproducibility and ecology can now be bridged thanks to the use of virtual reality (VR) for emotion elicitation. Headmounted displays (HMDs) provide a fully immersive experience by isolating the user from the external world and by inducing that immersive and engaging experience [13] lacking in IAPS, IADS or film-clips [14]. While tools such as IAPS allowed the processing of standardized emotional stimuli in different contexts (e.g., in hypoxia [15], the virtual experience allows for an investigation not only of discrete emotional stimuli, but also of emotional-inducing situations that contextualize them. Therefore, the use of VR in affective protocols has been shown to be an effective tool for basic and complex emotion elicitation [16], [17]. In this study, we propose 1) a novel, ecological fear induction VR scenario and 2) a rigorous computational model to estimate subjective feeling of fear as a function of specific components of skin conductance. To this aim, we have designed three immersive VR environments where the subjects could freely move into neutral, relaxing and fearful scenarios (preliminary results are shown in [18]). During their virtual experience, subjects' SC signals have been recorded using wearable sensors. To investigate fast-varying and slow-varying physiological correlates of fear, we have combined a robust and rigorous model for SC decomposition and feature extraction with a mixed-effects analysis at both single-and group-subject levels.

Experimental Protocol
We recruited 18 healthy volunteers (12 females) aged between 18 and 35, none of whom suffered from heart diseases or mental disorders (including phobias). Each participant gave the informed consent to take part in the study and was asked to report their interest in horror content, their experience level with VR systems, and to fill in the S.T.A.I. Y-1 questionnaire immediately before the experimental session to exclude participants with a high anxious state (i.e., scoring higher than 45). The bioethical committee of the University of Pisa (n. 14/2019) approved the experimental protocol.
We designed an experimental protocol consisting of four different sessions, each set in a different virtual reality scenario. To this aim, we adopted a commercially available portable head-mounted device: the Oculus Rift S (Lenovo Technologies and Facebook Technologies, USA). This headset presents a 2560x1440 LCD panel display, refresh rate up to 80 Hz, embedded tracking system with 6 degrees of freedom, and integrates headphones to reproduce sounds. The HMD was connected and supported by a PC to render and update the virtual contents. Throughout the whole experiment, each subject was equipped with a Shimmer3 GSR+ unit mounted on the wrist of the non-dominant hand. The Shimmer3 GSR+ unit is a wireless and wearable device able to record the SC simultaneously with photoplethysmography and accelerometer signals. More specifically, the SC was acquired at a sampling frequency of 50.33 Hz by placing two dry electrodes on the proximal phalanges of the index and medium finger of the non-dominant hand. To reduce the risk of artifact or spurious SC variations, each subject was asked to keep the non-dominant hand as steady as possible and not to speak throughout the experiment. Moreover, we limited any external stimulus to intensify the sense of presence.
The timeline of the experiment (see Fig. 1) was the following: Session 1: Relaxation R1. This session presented a blank space for 3 minutes and was designed to make the subjects relax and record their SC baseline level. Session 2: Neutral N. The subjects were immersed in a virtual living room for 3 minutes without any dynamic stimuli, to train and adapt them to the motion in virtual reality. Session 3: Fear F. During this session, the subjects visited for 5 minutes a scenario designed to elicit fear. To investigate the physiological reaction to fearful stimuli, we introduced some typical horror VR contents (see details in Section 2.2). Subjects were free to move within the scenario. According to the portion of the virtual environment explored, the subjects faced a certain number of such fearful content. Session 4: Relaxation R2. In the last session, we presented the same blank scenario as for R1 for 3 minutes. Since both the fear perception and the dynamics of the SC signals could be potentially affected by motion-sickness [19], we designed the virtual experience in order to limit the possible occurrence of such an issue (see Section 2.2). In addition, the subjects were asked to stand for the entire duration of the experiment, making the experience more natural, which is essential for reducing the motion-sickness risk. For the same purpose, the duration of the entire in-HMD session was set to less than 15 minutes, while the subject was free to explore the virtual scenario for just 8 minutes (i.e., N, F). At the end of each experiment, the subjects were asked to evaluate the perceived fear through a self-assessed fear score (SASF) consisting of a 10-points Likert scale ranging from 1 (not at all scared) to 10 (very scared). Using this scale, each participant assessed both the level of fear conveyed by the whole environment and the intensity of the fear felt for each discrete fear stimulus encountered within the scenario. In addition, we asked participants to report the possible onset of even light symptoms related to motion-sickness. All the eighteen participants did not report any symptom of motion-sickness.

Virtual Environments Design
The three virtual environments (see Fig. 2) were developed using the cross-platform game engine Unity3D (Unity3D 2019.4.8f1, Unity Technologies, USA). We selected Unity for both the flexibility in building virtual simulations and the ease of integration with the most popular VR HMDs. The 3D models representing the virtual contents in our experimental protocol were selected from the Unity Asset Store and free web services.
This immersive virtual setup is intended to maximize the effectiveness in inducing a psychological state comparable to what experienced in real settings. To this aim, a high engagement and sense of presence represents crucial aspects [13]. Accordingly, to improve the realism of the virtual setup, we designed the VR scenarios mimicking a first-person experience. This was obtained by using the locomotion system and the first-person vision provided by the Oculus Avatar Prefab, within the Oculus Integration package for Unity. Additionally, we further controlled the experimental conditions representing a key factor when dealing with emotional analyses and ANS (a-specific) signals. Accordingly, we limited the speed of participants' movement within the scenarios to avoid the insurgence of motion-sickness symptoms.
To provide active interactions between the virtual environment and users, we appropriately animated each of these VR contents through Mixamo, the Animator tool of Unity, and custom C# scripts.
To enhance emotigenic effects, the virtual environments were designed by taking into account several design parameters such as lighting, colors and geometries. Indeed, environmental light high temperature and intensity as well as direct lights may increase the arousal, with a proven effect also at neuro-physiological level [20]. Also, the tone of colors chosen in the VR scenes may influence the arousal level: warm colors tend to decrease arousal whereas cold ones have an opposite effect [21]. Finally, complex geometries determine a significantly higher arousal compared to the simple ones [22].
The fearful environment was designed according to scientific guidelines commonly applied to horror movies and games [23], aiming at eliciting the two main fear-inducing mechanisms shared among most individuals, i.e., fear of the unknown and jumpscares [24]. The effects induced by both visual and auditory cues have been broadly investigated. Most of the previous studies have focused on visual aspects (e.g., darkness, creatures belonging to the fantasy world) since fear reactions are mainly triggered by the presentation of visually perceptible fear-relevant stimuli [25], [26]. On the other hand, auditory cues have been also proven to heavily influence the emotional experience through suspenseful environmental sounds, voices and unknown audio [27]. According to the aforementioned guidelines, we designed the fearful environment to trigger both sensory channels. In particular, inspired also by commercial horror games, we set the fear environment in an abandoned hospital of 550 m 2 characterized by six rooms connected by corridors. Flickering dimmed lights were set to favor the darkness, and the related sound was included to generate suspenseful background audio. Without daylight and strong illumination sources, we allowed participants to light their paths and move through the darkness using only a virtual flashlight bound to the non-dominant hand. Finally, we included specific horror content that participants personally and actively faced at certain exploration points to stimulate punctual event-related fear responses. In particular, disfigured humans, zombie-like fantasy creatures, animals and sounds coming from doors opening and human voices were included in the VR environment given their demonstrated high fear elicitation capabilities [28]. To further confirm the fearful effect of the VR environment, the subjects were asked to score the induced fear level (i.e., SASF).
The two relaxation scenarios were designed according to the same aforementioned criteria. To minimize their emotional content, they were just blank scenes on a neutral, warm background color to avoid arousal stimulation. At the beginning of the R1 and R2 sessions, a black text suggested the subject relax and then it disappeared. To design the neutral scenario, we programmed a virtual living room representing neither a threatening nor a fearful situation. Unlike R1 and R2, during the N1 session, subjects could freely navigate the environment by using an Oculus touch controller. To make the virtual living room induce a neutral emotional effect, we dealt with light conditions and colours. Therefore, we combined both direct and indirect lights by setting the colour to a neutral temperature. Moreover, the intensity and range of visibility were set up to enable a vivid vision of the whole virtual room. Likewise, to choose the colours of the 3D objects that populated the living room we followed the same strategies applied for the colour of the lights.
Once designed, the three virtual environments were packaged in an application developed by combining the Unity platform and Visual Studio Code. The Unity application runs the VR scenarios at a fast frame rate of 78-80 Hz with a high rendering resolution of 3296 x 1776. The audio was reproduced in stereo mode at the maximum volume allowed by the integrated earphones. The choice of such system set-up parameters (e.g., frame rate, resolution, and audio volume) can be critical to induce that sense of immersion required to evoke realistic emotions (e.g., fear) in a simulated world [25]. In our system, the combination of high resolution and frame rate provided the participants with a fluid experience characterized by the absence of lags. Both these aspects are relevant to investigating a physiological response that is reliably elicited by a fearful experience [25]. Indeed, poor resolution can decrease the realism of the virtual environment by reducing the sense of immersion and distracting the subject from the task, thus reducing the fear elicitation effect [29]. On the other hand, the frame rate can be related to motion sickness due to the slow transitions between two consecutive frames [30]. Our experimental setting has successfully mitigated these risks. In addition, we adopted timed high-volume audio cues since they were demonstrated to be the best sound design to cause fear in subjects [27]. The application handled the experimental timeline, the fear contents spawn, and saved the event timestamps.

SC Processing and Decomposition
The SC signal consists of two different components characterized by different timescale and physiological mechanisms: the skin conductance level (SCL) and the skin conductance responses (SCRs). The SCL represents the slow-varying baseline of SC signal, and it is related to the general psychophysiological state of the subject [31]. The SCRs represent fast-varying superimposed SC changes. They can be spontaneous or can arise in response to an external stimulus with a latency of 3-5 seconds [32]. Such fast-varying fluctuations, which cannot be directly related to external stimulation, are called nonspecific SCRs (NSSCRs). Each SCR (including NSSCR) is the result of a neural burst in the sudomotor nerve activity (SMNA), which elicits the sweating from the glands under the skin surface of the hand. Therefore, both the components bring information about the autonomic nervous system functioning and are suitable to study the sympathetic response to fearful stimuli. From a psychophysiological point of view, SCL and SCRs represent correlates of emotional salience of situations and discrete stimuli respectively [33]. To decompose the SC signal into its components, we adopted the cvxEDA model [34], a decomposition algorithm based on Bayesian statistics and convex optimization. Together with the SCL and SCR signals, the model can estimate also the SMNA spike train signal. Before applying the cvxEDA algorithm, the SC signals were normalized using a Z-score transformation to facilitate model convergence.

Feature Extraction
To investigate the physiological responses to fearful stimuli, we conducted two separate analyses characterized by two different sets of features (Set1, Set2 in Table 1). Set1 was collected to identify the differences in the sympathetic state induced by the different VR environments, while Set2 was used to model the stimulus-specific physiological response to discrete stimuli (as those included in the fearful scenario) with different fear induction power.
Concerning Set1, we extracted the mean (SCLmean) and the standard deviation (SCLstd) of the SCL within nonoverlapped time windows of 20 seconds for all sessions. Moreover, we segmented the SCR component of each session with non-overlapped windows of 5 sec to analyze the NSSCRs. For each time segment we extracted the mean value of the NSSCR (NS-SCRmean); the number of NSSCR (NS-SCRn); the amplitude of the maximum NSSCR peak (NS-SCRPeak); and the sum of all peaks amplitude within the time windows (NSAmpsum). In addition, we included in the set the SC spectral power in the range 0.04Hz-0.15 Hz (EDAsymp) [35] computed within non-overlapped time windows of 60 seconds to provide the necessary frequency resolution. The extracted features within all-time segments were then averaged within the four experimental sessions (R1, N, F, R2).
Concerning Set2, we extracted features to characterize the stimulus-evoked response. Accordingly, we considered a 5-sec time window after the onset of each discrete fearful stimulus designed in session F, and we computed the mean (SCRmean), the maximum (SCRPeak) of the SMNA component; the number of SMNA peaks (SCRn), and the sum of their amplitudes (Ampsum); the latencies between the stimulus onset and the occurrence of the SMNA first peak (Lat1) and the SMNA with maximum amplitude (Lat2); the mean (SCLmean) and the standard deviation (SCLstd) of the SCL signal were calculated considering a longer 20-sec time window.

Statistical Analysis and Modeling
To evaluate significant SNS activity differences among the different experimental sessions, we conducted statistical analysis for each Set1 feature using nonparametric Friedman tests (p < 0:05) under the null hypothesis of no difference among the sessions. Post-hoc analysis between each pair of sessions was performed through a two-tailed Wilcoxon signed-rank test with Bonferroni correction. For each feature, we also tested possible significant interactions between subjects' gender and the reaction to the four experimental sessions. In particular, we applied a 2x4 mixed ANOVA design with gender as between factor (male, female) and experimental sessions as within factor (R1, N, F, R2).
Furthermore, we modeled the subjective fear quantified by the SASFs as a function of the SNS dynamics described by the SC features included in the Set2. More specifically, we performed an analysis based on linear mixed-effects (LME) models, where the SASFs was the dependent variable and the Set2 features were the independent fixed-effects. We adopted a mixed-effects model to account for the hierarchical structure given by the repeated measures across subjects and stimuli, thus exploring general linear relationships without violating independence assumptions of linear regressions. In addition, the mixed-effect architecture allows testing whether the differences among subjects (and, more interestingly, among fear contents) affects the reported fear perception.
First, the inter-class correlation coefficient (ICC) was computed to assess whether the variance of the SASFs depends on the hierarchical structure. Three different ICC were calculated: the SASFs variability explained by subjects, the SASFs variability explained by the discrete fearful stimuli, and the SASFs variability explained by the whole hierarchical structure. The statistical significance of the overall ICC measure was assessed by exploiting a simulation-based approach supported by a parametric bootstrap [36]. We used 1000 simulations to ensure a reliable proxy of the ICC confidence interval.
Second, possible collinearities between the Set2 regressors were tested to avoid instability in the model interpretation. In the case of collinear regressors, we selected the most correlated one with the SASFs. Once that such selection was performed, a group of models corresponding to the non-collinear SC-derived regressors in Set2 was further explored.
Third, to evaluate possible statistically significant effects of the regressors (i.e., the fixed-effects), we performed a multistep analysis: 1) As null model we used a random-intercept model without any fixed-effects to estimate the grouping variable effect on the SASFs intercept. This model served as a baseline model to test whether adding fixed-effects improved prediction reducing the error. 2) We built an LME model for each of selected Set2 features that were non-collinear. Specifically, each selected feature was added as a fixed-effect to the null model. LME models mathematically described the linear relationships between SASFs and each selected feature. In addition to the Set2 features, we tested the possible role of gender in explaining the SASFs. Accordingly, gender was included as a categorical fixed effect in the LME model. At this modeling stage, we also evaluated the addition of a multilevel variance structure for the fixed slopes, but we did not achieve a fitting convergence.
3) The contribution of each fixed-effect to explain the SASFs was statistically evaluated through pairwise comparisons between each LME model and the null model. We used the Akaike information criterion (AIC) as an index of the model quality. In addition, to test the statistical significance of each fixed-effect, we applied the likelihood ratio test (LRT), under the null hypothesis that there were no differences between each LME model and the null model. 4) All features with a significant LRT were included in a final LME model linearly combining multiple fixedeffects to explain the SASF variance. This full complexity model (FC-LME) describes the relationship between the subjective fear (SASFs score) and the selected SC-derived SNS activity. The statistical contribution of each fixed-effect and its interaction in the final FC-LME were analyzed using the same, aforementioned criteria. More specifically, we applied the LRT between the FC-LME and each LME model with a single fixed-effect that was significant in the step 3. Latency between stimulus and first SCR peak Lat2 Latency between stimulus and SCRpeak Fourth, a retrospective power analysis was accomplished to quantify the probability of the Type II error for each parameter of the FC-LME model. To this aim, we calculated the statistical power by relying on the effect size estimation (h 2 ). Hence, we computed the h 2 as the proportion of variance of the dependant variable associated with each fixed effect while controlling for the others.
Models fitting and comparisons were performed using the lme4 package of R [37].

RESULTS
In this section, we describe the results coming from both the statistical (Friedman) and the mixed-effects model analysis, applied to Set1 and Set2.
The answers to the questionnaires showed that for thirteen of the eighteen subjects the experiment represented their first VR experience, while five of them reported having an interest in horror content. All subjects got a score less than 45 on the S.T.A.I.-Y1 test, which meant a low/moderate level of anxiety state. Hence, none of them was excluded from the study.
On a scale from 1 to 10, the subjective fear related to the experiment as whole was rated 7.50 AE 1.50. On the other hand, discrete fear stimuli received a median score of 5.50 AE 1.50.

SC Correlates of Situational Fear
Friedman results on Set1 features showed significant differences in the SC dynamics (all p < 2:00e À 05). The following post-hoc Wilcoxon signed-rank tests performed for each feature revealed several significant variations between the pairs of experimental sessions, as shown in Fig. 3. Particularly, correlates of sessions F and R1 were significantly different for all the features. Likewise, when compared with R2, F session induced statistically differences in all the SC features except SCLmean. The comparisons between the N scenario and the other sessions were significant for most of the Wilcoxon tests. In fact, N session compared to R1 and R2 induced statistically significant differences in all the features with only two exceptions: SCLstd (N versus R1) and EDAsymp (N versus R2). Finally, when we compared the two relaxation stages (R1 versus R2), as expected, the significant differences decreased, and only SCL-related features (i.e., SCLmean and SCLstd) showed a significant p-value. This can be justified by the short timeline of the experimental protocol. In fact, the R2 stage can be considered a recovery session during which the SCL slowly return to baseline. In this case, our results showed a coherent significant decrease in the SCL activity from the F-related state to the R2-related state, but such a reduction trend did not have enough time to return to the baseline level (R1-related state) . It is worthwhile noting that more significant features were observed in the comparison between R1 and N than between N and F. This mostly concerns features related to the faster variations of the SC (e.g., NS-Ampsum, NS-SCRmean, NS-SCRn). This is due to the non-specific nature of the SCR considered in this Friedman analysis, which could be more affected by the transition from a state of immobility (R1) to one of virtual locomotion (N) than by the fearful punctual events that characterized the F session. This seemed to be confirmed by the remarkable differences that instead are shown by the SCLmean of N versus SCLmean of F. In fact, the SCLmean reflects changes in the general psychophysiological state of the subject (induced by F) more than an event-related response.
Moreover, the 2x4 mixed ANOVA analysis showed no significant interactions (p > 0:05) between the subject's gender and the differences in the SC features induced by the four VR experimental sessions.

ICC Analysis
The ICC estimating the subject-dependent variability and the stimulus-dependent variability were equal to 0.3681 and 0.182, respectively. Thus, the ICC of the variance explained by the whole hierarchical data structure (i.e., random-effects variance) was equal to 0:3681 þ 0:1820 ¼ 0:5501 with a 95% confidence interval ranging from 0.21 to 0.71 assessed by the non-parametric bootstrap analysis. Since the lower bound of the ICC confidence interval (i.e., 0.21) was greater than zero, we assumed both stimuli and subjects as random-effects of the LME models to capture significant fluctuations in psychometric fear perception due to such non-physiological factors.

Correlation Analysis
Significant collinearities were found between some of the Set2 features. Indeed, despite common in the SC analysis, Ampsum, SCRmean, SCRPeak and SCRn resulted significantly correlated among them (with a minimum correlation coefficient of 0.75). Since Ampsum showed the highest correlation with the dependent variable (i.e., SASFs), we included only this feature in the further LME analysis. Likewise, a significant correlation (r > 0:8) was observed between Lat1 and Lat2. According to the aforementioned strategy, only Lat2 was included in the model selection procedure due to its higher correlation with the SASFs. Thus, from the correlation analysis, the features considered for the LME modeling were: Ampsum, SCLmean, SCLstd, and Lat2.

LME Model With Single Fixed-Effect
For each Set2 feature selected after the correlation analysis, we fitted an LME model with a single fixed-effect (i.e., every selected feature). Table 2 shows the results of the comparisons, in terms of AIC and LRT, between every single fixedeffect model and the null model (i.e., the model with no fixed-effects).
In addition, we fitted a further LME model with the gender as single categorical fixed-effect. No significant relationship was found between SASFs and gender (x 2 ¼ 0:4282; p ¼ 0:5128).

FC-LME Model Definition
As final step of this analysis, Ampsum and SCLstd were linearly combined to build the FC-LME model. The same criteria of analysis used in the previous step were applied to compare the FC-LME model to each single fixed-effect models that was significant in the previous step. Table 3 describes the FC-LME model fitting and LRT results by reporting the slope value, the AIC, the x 2 resulting from the LRTs, and the corresponding p-value.
Both fixed-effects of the model, i.e., Ampsum (x 2 ¼ 9:32; p ¼ 0:0022) and SCLstd (x 2 ¼ 5:34; p ¼ 0:0200), were significantly associated with the SASFs outcome. The conditional r 2 , which indicates the variance explained by the model's fixed and random-effects, was equal to 0.53. The mean square error (MSE) was equal to 0.81. The marginal r 2 , associated with the only fixed-effects' capability to explain the SASFs variance, was equal to 0.30.
Finally, we evaluated the interaction term between the two selected fixed-effects (Ampsum and SCLstd). The FC-LME model including the interaction term was compared through LRT with the FC-LME model. This test revealed that the interaction term significantly contributes to the SASF prediction (x 2 ð1Þ ¼ 8:06; p ¼ 0:004) with a slope of -0.032. This final model showed a conditional r 2 equal to 0.597, an MSE equal to 0.83 and a marginal r 2 equal to 0.40.
The (1) formally describes the resulting FC-LME model Where SASF ij represents the j th self-assessed fear score of the i th subject. b 1 (0.016), b 2 (10.183) and b 3 (-0.032) are the fixed slopes corresponding to Ampsum, SCLstd and their interaction, respectively. m denote the intercept (2.952) of the model. a ð1Þ indicates the Gaussian distributed subject random-effect, and a ð2Þ is the Gaussian distributed stimulus random-effect. Focusing on the estimated model parameters, we obtained an effect size of 0.21 (CI=[0.08 -1]), 0.23 (CI= [0.09-1]), and 0.16 (CI=[0.05-1]) for Ampsum, SCLstd and the interaction between them respectively. In addition, based on the computed effect size, we estimated a statistical power equal to 0.95 for the two regressors and 0.87 for the interaction term. The model fitting results are also shown in Fig. 4 through a 2-D representation (on the right side). The Fig. 4 shows how the value of SASFs increases proportionally to both SCLstd and Ampsum and decreases as a function of the interaction.

DISCUSSIONS
In this work, we combined VR techniques and statistical modelling to explore the possible linear relationships between the intensity of subjective fear and the slow-and fast-varying autonomic response to fearful situations and stimuli. We presented a novel virtual environment for an ecological but controlled fear elicitation that was designed to combine strategies commonly applied to horror movies and games with techniques to successfully reduce the probability of motion-sickness insurgence. Our study confirmed the effectiveness of VR technology in inducing emotional states and in reducing the gap between a realworld emotional elicitation and highly controlled laboratory conditions. In fact, both the self-assessed level of fear and the analysis of the SC-derived features varied coherently with the different emotional salience attributed to each experimental session. Particularly, we observed an increment in the sympathetic neural burst frequency and intensity during the immersion in the VR fearful scenario. Moreover, considering the SCL-derived features, which allow inferences about the subject's slow-varying emotional processes, a significant variation of the SCLmean among the experimental sessions suggested an alteration in the subject's state towards a more aroused condition. The increase of the SNS activity was further confirmed by the statistical significance of the sympathetic activity index, EDAsymp, showing a higher activity during fear-inducing (F) sessions compared to relaxation (R1, R2) ones. Overall, the trend of the SC features among the VR sessions (see Fig. 3) indicated an alteration of the general subject's state toward a more aroused condition during the F session as confirmed by SASFs scores. This association between selfreported levels of fear and their physiological correlates corroborates the reliability of the proposed VR environment in eliciting the emotional state of fear.
These VR scenarios were used to develop a new approach to investigate the relationship between the subjective feeling of fear and the sympathetic nervous system activity, by analyzing of one of the most widely used ANS correlate (i.e., the skin conductance). To this aim, we proposed an LME approach to fully describe the hierarchical structure of the data given by repeated measures across subjects and stimuli. As expected the subjective perception of fear resulted significantly influenced by the common intrinsic differences among subjects, as well as by the different fearful contents presented in the virtual scenario. On the other hand, gender difference did not seem to significantly relate with the self-reported fear level. Quantifying the relationship between self-reported fluctuations in fear and ANS correlates (e.g., SC features) represents an essential step towards a more objective evaluation of emotional states, which can fail to be effectively recognized through introspection due to clinical (e.g., alexithymia [38]) or experimental (e.g., perceptually-or emotionally-subliminal [39] stimulation) conditions. The model design and selection procedures brought relevant insights concerning which SC features quantified the sympathetic dynamics involved in the processing of fearful scenes. In particular, a significant linear relationship was found between the subjective scores of fear perception and the sum of the sympathetic tone variation and the neural burst amplitudes induced by discrete fearful stimuli. Particularly, a variation in burst frequency and intensity (i.e., Ampsum) explained a significant variation in the reported level of subjective fear. This suggests that, among the SC features considered, these represent the most reliable correlates of self-reported fear in response to VR stimuli. The temporal properties of such correlates are coherent with the hypothesis of a quick-and-dirty, short path  coordinated with a longer, more precise path underlying the processing of fearful stimuli [1]. Indeed, the significant contribution made by SCLstd to the model could reasonably indicate the oscillation in the sympathetic tone induced by a wider difference in the perceived salience of stimuli. Thus, fearful stimuli could induce a twofold reaction that occurs on different time scales and involves different phenomena: a fast-varying response, represented by Ampsum, occurs more quickly as the direct manifestation of the fight-or-flight response; whereas a slow-varying response, represented by SCLstd, represents the variability in the perceived emotional salience of each VR scenario. This pattern fits with bifactorial models of fear regulation describing the relation between a faster (and not necessarily conscious) threat-detection circuit and a slower circuit allowing the cognitive emotional representation [1]. The different ability of these two factors to affect each other was proposed as a marker of disorders such as specific phobias [3], [39] or Panic disorder [40]. Indeed, both conditions are characterized by the inability to cognitively contain an exaggerated reaction of fear induced respectively by discrete stimuli and situations, even if the patients recognize such emotion to be irrational and unmotivated. Beyond this psychophysiological interpretation, the statistical significance of SCLstd could also reflect mathematical aspects related to the cvxEDA optimal solution of the quadratic program. More specifically, the decomposition process could attribute a part of the SCR-related variance of the SC signal to the SCL component. In this view, the LME model could have combined the split contributions which are instead the same representations of the faster SCR response.
As a final result of the study, selected features were linearly combined to build the so-called FC-LME model. Differently from commonly used machine learning and linear regression approaches, the LME context affords a flexible and straightforward solution for modelling highly structured data characterized by hierarchical clusters of non-independent observational units. Indeed, on the one hand, LME analysis allowed to overcome the limitations of learning models and generalized linear models without random-effects, which cannot analyze repeated measures across subjects and cluster-related variability without manipulating data to aggregate non-independent observations. On the other hand, a mixedeffects model's structure mitigates possible limitations due to the reduced number of participants (i.e., 18). However, modelling the idiosyncrasies in the data due to the randomeffects structure leads to a robust inference of the fixed-effects without inflating the Type II error rate. [41]. This was partially corroborated by the calculated statistical powers that demonstrated a high probability of having detected real effects of Ampsum and SCLstd on the subjective fear perception.
The FC-LME model could help to overcome the issues characterizing the standard self-assessed measures of emotion, including fear. Indeed, tools based on self-assessment measures are limited by subjects' hesitation in giving honest answers when it is socially undesirable or by the difficulty in labelling the emotions, introducing subject-dependent biases in the measures. Instead, the proposed FC-LME model provides a random-effects structure that explicitly accounts for the processes that generate such subject-depending variability. It hence represents an objective tool to integrate (or replace, where necessary) the self-assessed measures.
During the different processing and modelling phases we faced and addressed some potential critical issues. First, the discrete ordinal nature of the dependent variable (i.e., the Likert scale ranging from 1 to 10) could be critical due to the linearity of the model, which assumes Gaussian distribution of residuals and linear variation of the regressors across the ordinal levels of the dependent variable. However, the high number of levels (i.e., 10) of the ordinal variable mitigates the negative effect on the model fitting. Indeed, previous studies have confirmed the LME robustness with such a Likert scale [42]. Moreover, on the definition of the random-effects structure, we did not consider random slopes since the outputs of the LRT revealed no significant contributions to the variance explained by the model. However, despite being statistically plausible, it is unlikely that the slopes associated with biological regressors do not vary across subjects and stimuli. On the one hand, the limited number of observations could affect the possibility of successfully fitting a covariance structure that considers both random intercepts and random slopes across subjects and stimuli. On the other hand, the subjects were free to explore the fear scenario and could miss some fearful events, causing an unbalanced number of observations (i.e., discrete fear stimuli responses) for each subject. This might cause an unstable model with variance components collapsing to zero, especially for random slopes model fitting [43]. In the light of the above, although considering random slopes might improve the fixed-effects inference and control the Type I error rate [44], fitting a full-complexity covariance structure is often unrealistic for most biological datasets due to the volume of data required to estimate the random intercepts and slopes variances, and their covariance [45]. Accordingly, the proposed FC-LME model defines the most complex random-effects structure allowed by the available data [45]. Furthermore, a random intercept model that mathematically describes the idiosyncrasies within the data was less critical in terms of model convergence, once that the dependent variable levels are greater than five [46], [47].

CONCLUSION
The present study provides an innovative tool that uses VR to investigate both slow-and fast-varying autonomic correlates of fear in an environment combining a high realism with infinite possibilities of experimental manipulation. Analyzing the skin conductance features of healthy participants, we found a correlation with self-reported levels of subjective fear, thus allowing a more objective and direct access to this crucial information. If the association between skin conductance and self-reported levels of subjective fear will be confirmed and validated in clinical samples, it could represent an additional piece of information to be considered while diagnosing mental disorder characterized by pathological fear (e.g., specific phobias; PTSD; panic disorder). The same approach could be useful also to achieve a more objective assessment of symptom's severity, since psychophysiological correlates of fear could be less affected than self-reports by psychological variables. Furthermore, the possible dissociation between physiological and subjective outcomes could allow the characterization of complex clinical situations such as psychiatric comorbidities. In this regard, bringing specific phobias as an example: 1) an overlapping narcissist trait could induce to exaggerate symptoms severity up to pathological thresholds, but observing a little electrodermal response after the exposure to the phobic stimulus could help therapists to achieve a more accurate diagnosis; 2) an overlapping major depressive disorder could lead to exaggerate own powerlessness toward the feared stimulus, thus shifting phobic symptomatology up to clinical levels; 3) conversely, an overlapping manic phase of bipolar disorder could induce to underestimate the severity of own symptoms. Additional future studies will focus on corroborating the relationships between the levels of selfreported fear and the ANS correlates. To this aim, we will go towards a continuous scale of the SASFs increasing also the number of discrete fearful stimuli within the VR environment, and the number of participants. Furthermore, we will also explore the influence of further demographic data (e.g., age, sex, interest in horror content) on fear perception by improving the subjects' profiling and including such information as covariates in the model. Beyond this further model validation, a multimodal, multivariate, multiorgan and multiscale [48] LME model will be evaluated by integrating physiological features from different ANS correlates such as the heart rate variability and the respiratory signals.