Are you afraid? Neurophysiological responses to multidimensional fear in an interactive semi-immersive environment

1 Introduction

Fear recognition is critical for both clinical evaluations, such as trauma assessment, and for understanding normative emotional processing (Ayuso-Moreno et al., 2020; Ribas et al., 2018; Russo et al., 2024; Battaglia et al., 2022). As a fundamental defensive state, fear mobilizes the central and peripheral nervous systems to facilitate rapid detection and response to threats, a reaction commonly known as the fight-or-flight mechanism (Cannon, 1915). This state manifests through increased autonomic arousal and altered cortical activity, whose physiological measurement enables objective assessment and distinction between adaptive vigilance and pathological conditions such as anxiety isorders.

Electroencephalography (EEG), in particular, has been widely used to study rapidly changing cognitive processes with high temporal resolution (Berger, 1929; Luck, 2014), including mental fatigue, vocational interest, cognitive load, attention, and emotions (Ramírez-Moreno et al., 2021; Candela-Leal et al., 2025a,b, 2026; Blanco-Ríos et al., 2024). This is primarly because the spectral EEG features are associated with underlying cognitive states, such as: Delta (δ, 1–4 Hz) to deep sleep, theta (θ, 4–8 Hz) to working memory, alpha (α, 8–12 Hz) to attention, beta (β, 13–30 Hz) to executive action, and gamma (γ, 30–100 Hz) to higher-order cognition (Kropotov, 2009). In the current study, only the theta, alpha, and beta bands will be analyzed, given their roles in mental effort (Candela-Leal et al., 2025b; Villafaina et al., 2021b), attention (Benedek et al., 2014; Ray and Cole, 1985; Klimesch et al., 1999; Jensen, 2002; Sauseng et al., 2005; Rihs et al., 2007), and executive action (Lamarche and Ogilvie, 1997; Pineda-Hernández, 2022), respectively.

Complementary to brain measurements, autonomic activity can be assessed through signals regulated by the Autonomic Nervous System (ANS). In particular, heart rate and sweat gland activity enable the extraction of physiological features such as Heart Rate Variability (HRV) and Electrodermal Activity (EDA), which reflect sympathetic and parasympathetic modulation (Wehrwein et al., 2016; Ghiasi et al., 2020). HRV is based on the successive heartbeats variation over a time interval (Thayer et al., 2012), in which the cardiovascular center adjusts the heart rate by shifting between sympathetic and parasympathetic nervous systems (Loewi, 1921; Levy, 1971). A lower level of HRV, sympathetic activation, is related to reduced regulatory capacity to adapt to challenges or cognitive stressors (Shaffer et al., 2014), while higher HRV is an indicator of resilience, behavior flexibility or cognitive load (Beachaine, 2001; Berntson et al., 2008). EDA measures changes in the electrical conductance of the skin resulting from the modulation of sweat gland activity by the sudomotor nerves (Ghiasi et al., 2020); these changes are determined exclusively by the sympathetic nervous system. Therefore, EDA can be used to reliably measure sympathetic function, as well as cognitive arousal (Ghiasi et al., 2020; Posada-Quintero and Chon, 2020). The most salient characteristic of EDA signals is the appearance of Skin Conductance Responses (SCRs), which are phasic components from a sympathetic reaction to stimuli. Besides the analysis of rapid and smooth changes trough SCRs, the Skin Conductance Level (SCL) measures slow shifts in the signal, defined as the overall conductance obtained from the tonic component of EDA. In conjunction, SCRs and SCL assess the subject's response in event-related experiments or tonic stimuli tests, such as changes in condition, workload, or cognitive stress (Posada-Quintero and Chon, 2020; Posada-Quintero and Bolkhovsky, 2019).

At the brain level, human functional Magnetic Resonance Imaging (fMRI) studies have provided insights into neural circuitry related to fear, including the amygdala, hippocampus, and prefrontal cortex (Åhs et al., 2015; Picó-Pérez et al., 2019; Milad and Quirk, 2012). On the other hand, EEG studies have identified frontal (Rezaei et al., 2025; Kim et al., 2019; DeLaRosa et al., 2014; Apicella et al., 2023; Demaree et al., 2005), parietal (Rezaei et al., 2025; Ugalde-Canitrot et al., 2025; Alemayehu et al., 1995; Yaşgüçlükal et al., 2021; Okazaki et al., 2008), temporal (Ribas et al., 2018; Byeon et al., 2023), and occipital (DeLaRosa et al., 2014; Ugalde-Canitrot et al., 2025; Okazaki et al., 2008) lobes as cortical structures related to fear; relating theta (Rezaei et al., 2025; Kim et al., 2019; DeLaRosa et al., 2014), alpha (Russo et al., 2024; Kim et al., 2019), and beta (Ribas et al., 2018; Russo et al., 2024; Rezaei et al., 2025; Kim et al., 2019; Apicella et al., 2023; Okazaki et al., 2008; Byeon et al., 2023; Bălan et al., 2020a) frequency bands. These have different physiological meanings depending on the oscillatory frequency and brain area. For example, in response to threatening images, there is a theta burst in the occipital region, which processes visual information, followed by a theta burst at the frontal region for decision-making and planning (DeLaRosa et al., 2014). Conversely, beta activity in the same region increases in situations that require significant attention and in which a perceptual shift occurs (Rezaei et al., 2025; Kim et al., 2019; Okazaki et al., 2008).

Beyond brain dynamics, autonomic responses have also been used as objective markers of fear, mainly in the study of fear conditioning. In this context, sympathetic arousal is commonly inferred from SCRs (Staib et al., 2015; Castegnetti et al., 2016), i.e., threatening or fear-relevant cues elicit larger and more frequent SCRs than neutral or safe cues in Pavlovian fear conditioning and extinction paradigms (Christopoulos et al., 2019; Marin et al., 2017). However, SCR analysis is limited by its inability to capture parasympatheti, vagal, activity. This limitation is relevant, as fear-conditioned stimuli have been shown to elicit transient heart rate deceleration, a response generally believed to reflect predominantly vagal activation (Castegnetti et al., 2016). Consequently, increasing efforts have focused on the use of HRV as an index of fear conditioning, given its potential to assess sympathetic and parasympathetic activation in response to stimuli (Battaglia et al., 2022; Battaglia and Thayer, 2022).

While traditional neurophysiological studies of fear often rely on highly controlled, artificial laboratory settings (Ghiasi et al., 2020; Thayer et al., 2012; Staib et al., 2015; Castegnetti et al., 2016), immersive and semi-immersive technologies offer more ecologically valid contexts for studying emotion, with semi-immersive environments enabling free expression and movement while also enhancing immersion (Kisker et al., 2021). The effectiveness of Virtual Reality (VR) in eliciting fear is highly modulated by presence (Riva et al., 2015), which is the psychological sense of “being there” in the Virtual Environment (VE) (Witmer and Singer, 1998), characterized by spatial presence, involvement, and realness (Schubert, 2003). Presence is enhanced by immersion (Diemer et al., 2015), which depends on the equipment's characteristics used to create and display the VR (Gromer et al., 2018), such as display quality and real-time feedback (Slater, 1999), and has been associated with increased frontal beta activity (Tychkov et al., 2023). Therefore, multi-modal interfaces are essential to provide multisensory stimulation, improving fidelity and presence (Cooper et al., 2018; Honegger et al., 2021; Riva et al., 2007); without them, unimodal interfaces often feel hollow and unrealistic (Maggioni et al., 2019; Mishra et al., 2020).

Consequently, given that VR can evoke emotions (Higuera-Trujillo et al., 2017; Chirico and Gaggioli, 2019; Kisker et al., 2021), they are increasingly used in closed-loop phobia therapies. These systems monitor physiological signals to dynamically adapt exposure scenarios in real-time (Bălan et al., 2020a,b). While exposure may acutely induces stress (EDA and HR) during the experience (Moldoveanu et al., 2023), it reduces post-treatment anxiety across various phobias (Bălan et al., 2020b; Diemer et al., 2016), including acrophobia (fear of heights) (Russo et al., 2024; Apicella et al., 2023; Bălan et al., 2020a; Freeman et al., 2018), ophiophobia (fear of snakes) (Toma et al., 2020; Go et al., 2020; Ding et al., 2020), claustrophobia (fear of enclosed spaces) (Rezaei et al., 2025), pyrophobia (fear of fire) (Nazarova et al., 2023; Obukhov et al., 2023) and fear of public speaking therapy (Reeves et al., 2021; Seuling et al., 2024).

The current paper addresses the static settings and unimodal limitations by integrating non-invasive multi-modal neurophysiological (EEG, BVP, EDA) recordings of participants at a real-time interactive (psychophysiological-based) immersive environment, considering an 8-electrode EEG (OpenBCI) at theta, alpha, and beta spectral bands; in addition to a biometric bracelet (Empatica E4) to dynamically modulate the surrounding projection-based VR setting, including visualizations, lighting, and sound. To identify key neurophysiological responses related to fear-based environments, a pairwise analysis with baseline and an interscene analysis with other expressive scenes within the immersive environment is considered. Specifically, addressing the research questions:

• How do neural and autonomic responses differ between baseline and fear-related scenes in a semi-immersive environment?
• How do neural and autonomic responses change across different semi-immersive interactive fear-related scenes?

The hypothesis is that participants in the immersive fear scenario will exhibit increased stress-related autonomic responses (EDA, HRV). On the other hand, in all immersive scenes, there should be increased brain activity associated with visual information processing (occipital), planning (central), emotional processing (parietal), and decision-making (frontal). Thus, showing that fear-related semi-immersive scenarios elicit information processing and activate the sympathetic nervous system at both psychological and physiological levels.

The paper is organized as follows: Section 2 describes the participants (Section 2.1), the signal acquisition devices and process (Section 2.2), the experimental protocol followed (Section 2.3), in addition to also detailing the EEG, BVP, and EDA signal pre-processing (Section 2.4) and processing (Section 2.5), as well as the statistical analysis (Section 2.6). Section 3 uses the normalized EEG, HRV, and EDA features to assess shared brain and autonomic activation patterns across scenes (Section 3.1), and to show cross-device relationships (Section 3.2). The most important findings are discussed in Section 4, where Section 4.1 addresses research question (1) and Section 4.2 addresses research question (2). In addition, Section 4.3 adds complementary inter-device findings on both research questions. Finally, Section 4.4 discusses the study's limitations and Section 5 concludes with the relevance of this article.

2 Materials and methods

The following sections describe the participants, experimental procedures, and analytical methods used to evaluate the neurophysiological responses associated with the fear experience.

2.1 Participants

Participants consisted in 10 healthy adults (M = 4, F = 6), with ages ranging from 18 to 36 years (μ = 25.11, σ = 6.21). Participants signed an informed consent form in accordance with the 1964 Helsinki Declaration. They were informed about the study, their rights, risks, and benefits, and remained in close contact with the research team. The exclusion criteria ruled out participants who were undergoing pharmacological treatment or had cardiovascular or physical conditions that limited their ability to perform the required movements. The experimental research protocol was granted an exemption (protocol number EHE-2023-03) by the Institutional Research Ethics Committee (CIEI) of Tecnologico de Monterrey. The experiment and signal acquisition were made in collaboration with Centro de Estudios Críticos en Cultura Contemporánea (CECRIT), Universidad Autónoma de Queretaro (UAQ). The space measured 4 by 4 m, and data were recorded over 1 week (29-Apr-2024 - 03-May-2024).

2.2 Signal acquisition

EEG signals were acquired through the Ultracortex Mark IV (OpenBCI, New York, NY, USA), a dry-electrode reconfigurable headset that includes eight channels: Two at the prefrontal region (Fp1, Fp2), two at the central region (C3, C4), two at the parietal region (P7, P8), and two at the occipital region (O1, O2); in addition to two ear clips on both earlobes. EEG data was measured in microvolts (μV) at a sampling frequency of 250 Hz. OpenBCI GUI (OpenBCI, New York, NY, USA) was used to measure the impedance of electrodes and electrodes were screwed more tight to the participant's scalp until the impedance values on most channels were less than 750 kΩ (Romo-De León et al., 2024), indicated by a green electrode color at the interface, this software was also used to transmit EEG scalp data from the Cyton board to the computer through Bluetooth by leveraging the OpenBCI dongle.

Autonomic responses were acquired through the Empatica E4 wristband (Empatica, Cambridge, MA, USA), a wearable device designed for ambulatory physiological monitoring. The E4 is capable of recording EDA at a sampling frequency of 4 Hz, skin temperature at 4 Hz, triaxial acceleration at 32 Hz and Blood Volume Pulse (BVP) at 64 Hz; thereby allowing the derivation of Interbeat Intervals (IBIs) and Heart Rate (HR). Data were transmitted wirelessly in real-time to the host computer via a second dedicated Bluetooth dongle, enabling continuous data acquisition during the experimental protocol.

2.3 Experimental protocol

The experimental design (Figure 1) consisted of a learning experience divided into four blocks, referred to as “scenes,” each consisting of specific 3-min tasks. Participants undertook each scene sequentially, starting with a 2-min eyes open baseline resting period, then scene 1, 2 and 3, and ending with scene 4. Each scene was separated by a 1-min washout period, during which the investigator explained the upcoming task, saved the previous recording, and initiated the next. Scenes order was fixed given the constructive nature of the experience, where decisions in earlier scenes (e.g., words chosen in Scene 2) affect subsequent scenes (e.g., words spoken in Scene 3). The scenes' design was based on (Romo-De León et al., 2024, 2026), which is grounded in Kant's philosophy of the Theory of Knowledge and follows a structure of cognitive learning and integrative levels. Each scene had an underlying logic that provided it with a fundamental structure, representing different levels of cognition into which modern philosophy has divided experience, according to Kant, the primary reference in the modern theory of knowledge (Dennis, 2020; Kuplen, 2023; Webb, 2020).

More particularly, scene 1 (Embodied Fear) consisted of the physical experience, in which no language is involved, only bodily expression. Participants were exposed to a stressful situation, evoked through images and sounds associated with the experience of fear, which were iteratively co-created by creative artists and the research team to build an immersive environment for closed-loop neurofeedback experience, rather than using standardized affective stimuli. Participants were invited to react physically to the environment, while their movements altered the distance of the projection and generated specific sounds as they moved their limbs, thereby creating a semi-immersive experience. Throughout all scenes, EEG and autonomic signals were recorded to form the dataset analyzed offline in this study, utilizing Unix timestamps for cross-device synchronization. Simultaneously, a parallel stream of these signals was processed online using the pipeline described by Ortiz-Espinoza et al. (2025) to drive the adaptive audiovisuals. At the same time, an algorithm calculated a real-time categorical fear score based on an EEG-based Valence-Arousal-Dominance (VAD) emotion model (Blanco-Ríos et al., 2024), previously trained on the Database for Emotion Analysis using Physiological Signals (DEAP) (Koelstra et al., 2012). Following Bălan et al. (2019), VAD values were projected into a 3D emotional space to determine four fear levels across the hate-overconfident diagonal: (1) no fear [V: (7, 9), A: (1, 3), D: (7, 9)], (2) low fear [V: (5, 7), A: (3, 5), D: (5, 7)], (3) medium fear [V: (3, 5), A: (5, 7), D: (3, 5)], and (4) high fear [V: (1, 3), A: (7, 9), D: (1, 3)]. This score was displayed to participants in real-time and modified the intensity of the environment across scenes, making it more intense or more relaxed depending on whether a higher or lower fear level was detected, respectively. Similarly, the flashing of LED lights accompanied these detected ranges by blinking with greater or lesser intensity.

Scene 2 (Conceptual Fear) consisted of a reflexive exercise in which participants presented a conceptual explanation of their previous experience. On the screen, a cloud of words associated with the emotion of fear was projected. The participants said aloud into a microphone the words that best represented the experience of fear, or proposed alternative words. These suggestions were generated using a recommendation system based on word representations extracted from a predefined corpus. The corpus treats the emotion of fear, comprising contemporary novels, newspaper articles, social media posts, and essays.

Scene 3 (Discursive Fear) involved applying the previous concepts in articulated sentences, going beyond the conceptual sphere of knowledge into the discursive one, where nuances such as meaning and emphasis were taken into account. The screen displayed the text fragments that most closely matched the words selected in Scene 2, based on Jaccard similarity. Each paragraph consisted of two or three sentences. The participant was then asked to read the fragments out loud into the microphone in the most emphatic and comprehensible manner.

Scene 4 (Integrative Fear) ends with a symbolic synthesis of the whole experience, appealing to figurative images that can embody different meanings and interpretations of the same fact. The participants were seated, and a camera recognized their facial features through MediaPipe. The Vanitas painting was projected, and they observed it and explored the emotions it evoked by engaging facially. This was done online through GenAI by leveraging a NVIDIA GeForce RTX 3090 eGPU (24 GB VRAM), explained in greater detail in (Ortiz-Espinoza et al., 2025).

2.4 Signal pre-processing

EEG data were pre-processed by first assigning the respective EEG channel positions (Fp1, Fp2, C3, C4, P7, P8, O1, O2). Performing the standarized pre-processing (PREP) pipeline (v0.57) (Bigdely-Shamlo et al., 2015) via MATLAB R2024a (MathWorks, Natick, MA) and EEGLAB toolbox (v2024.0) (Delorme and Makeig, 2004), which detrends using a high-pass filter (1 Hz), removes line noise by using multiples of 60 Hz, and re-references the EEG data through a robust average reference, which iteratively identifies bad channels in order to interpolate them using the other channels. Furthermore, a Finite Impulse Response (FIR; v2.8; 4–30 Hz) was applied; later, Artifact Subspace Reconstruction (ASR) algorithm (v2.0) (Mullen et al., 2013) (κ = 20) was implemented, which improves signal quality by correcting movement-related artifacts (Cruz-Garza et al., 2020); finally, an Independent Component Analysis (ICA) was performed to classify and remove the non-brain components, including line noise, eye, heart, and muscle, with ICLabel (v1.6).

Using this EEG pipeline, all windows were retained, as ASR and PREP automatically managed artifact correction and bad channel interpolation, respectively. Non-brain components were identified spatially, spectrally, and temporally using ICLabel; those with ≥70% brain probability were retained, removing a maximum of two non-brain components per participant. Analysis was restricted to complete datasets (calibration plus four scenes). Out of 10 subjects, two were excluded: one due to a corrupted calibration file (S5), and one due to rank-deficiency of 1 during ICA caused by channel interpolation (S1); with an effective sample of eight subjects.

BVP signals were denoised using an adaptive frequency-domain approach based on partial Discrete Fourier Transform (DFT) reconstruction, as presented in (Sadhukhan et al., 2018). For each recording, the raw BVP signal was segmented into 10 s windows with 2% overlap. Low- and high-cut frequencies were estimated from the first signal segment. The low-cut frequency was derived from HR estimation with HeartPy toolkit (v1.2.7) (van Gent et al., 2018) following high frequency noise attenuation with an 8 Hz Butterworth low-pass filter. The high-cut frequency was determined from the signal spectrum as the highest frequency whose magnitude exceeded 15% of the dominant spectral peak above the low-cut threshold. Frequencies outside of the estimated band of interest were then rejected from each segment, and signal reconstruction was performed via inverse DFT. Finally, the preprocessed BVP signal was obtained by concatenating the reconstructed segments, discarding overlapping sections.

EDA signals were preprocessed using a data-driven empirical iterative method adapted from (Gautam et al., 2018). For each recording, the raw signal was iteratively smoothed. At each iteration, local maxima and minima were detected and midpoints between paired extrema were computed. A cubic spline was then fitted through these midpoint values to obtain a smoothed component, which would then undergo the same treatment. This procedure was repeated for up to seven iterations, as suggested in the original methodology, or stopped earlier if insufficient extrema pairs were detected. Finally, the last smoothed component would be substracted from the original EDA signal to obtain a residual to be used for subsequent analysis.

2.5 Signal processing

Pre-processed EEG brain signal was processed employing Python (v3.9.0) and Brainflow (v5.6.0). Using a frequency analysis based on the Fast Fourier Transform (FFT), a 1-s PSD with a half-second overlap was calculated. Afterwards, the average power for each of the following band ranges was calculated: Theta (θ, 4–8 Hz), alpha (α, 8–12 Hz), and beta (β, 12–30 Hz). Afterwards, the averaged pooled brain activity was calculated for each brain region: Frontal (FP1, FP2), Central (C3, C4), Parietal (P7, P8), and Occipital (O1, O2).

In addition, taking into account the combined average brain activity by region, frequency band ratios were also calculated, considering an average of theta at frontal (Fp1, Fp2) and central region (C3, C4), an average of alpha at parietal region (P7, P8), and an average of beta at frontal region (Fp1, Fp2), according to Kamzanova et al. (2011). The EEG-based indices calculated were: Engagement index (Equation 1a), considering cognitive processing in contrast to a more passive state (Pope et al., 1995); Excitement index (Equation 1b), defined as high alert and attentive (Pfurtscheller and Lopes da Silva, 1999); and Fatigue Index (Equation 1c), with mental weariness and required at sustained attention (Ramírez-Moreno et al., 2021).

Engagement(θ,α,β)=βθ+α(1a)

Excitement(α,β)=β/α(1b)

Fatigue(θ,α)=θ/α(1c)

Pre-processed Empatica E4 data (BVP and EDA) were segmented into 90 s windows with 50% overlap. Feature extraction was performed independently for each participant and condition, and features were subsequently concatenated across modalities within aligned windows. Given the limited length of the available recordings, HRV features were obtained using an ultra-short-term analysis approach. Accordingly, a window length of 90 s was selected based on recommendations summarized in Shaffer and Ginsberg (2017). This duration is primarily suitable for estimating metrics that reflect rapid beat-to-beat variability, particularly the root mean square of successive RR interval differences (RMSSD) and HF power. Nevertheless, HRV estimates derived from 90 s segments should not be regarded as interchangeable with those obtained from conventional 5 min recordings, particularly for metrics that require longer observation periods (Shaffer and Ginsberg, 2017; Baek et al., 2015; Munoz et al., 2015; Shaffer et al., 2020). To maintain temporal consistency across modalities and facilitate window-level feature integration, EDA signals were segmented using the same 90 s window length. A 50% overlap produced a 45 s step, allowing the initial, intermediate, and final portions of each condition to be represented while balancing temporal resolution and redundancy. Although no physiologically established overlap percentage exists for HRV or EDA, overlapping windows have been used in continuous cardiovascular and electrodermal monitoring (Hovsepian et al., 2015; Gjoreski et al., 2016; Schmidt et al., 2018; Nath and Thapliyal, 2021). Using the aforementioned pipeline, all windows were retained for analysis. Out of 10 subjects, two were excluded due to corrupted files (S7, S10); with an effective sample of eight subjects.

For BVP, peak detection and inter-beat interval (IBI) series estimation were conducted using the HeartPy toolkit (van Gent et al., 2018). For this, a smoothed version of the signal (resulting from the application of a moving average filter) was used as a reference for peak detection. Detected peaks were then verified and corrected, and the resulting RR intervals were cleaned from outliers using the quotient filter method (Piskorski and Guzik, 2007). From the corrected RR series, heart rate (to be referred to as BPM from here onward) as well as HRV features were obtained for each time window. The computed time-domain HRV metrics included the standard deviation of NN intervals (SDNN), root mean square of successive RR interval differences (RMSSD), and the percentage of successive RR intervals that differed by more than 50 ms (pNN50). Frequency-domain metrics were estimated from the power spectral density of the corrected RR series using Welch's method, with LF power (0.04–0.15 Hz), HF power (0.15–0.40 Hz), and total power computed by numerical integration of the corresponding spectral bands. The LF/HF ratio was subsequently calculated from the estimated LF and HF powers (Shaffer et al., 2014). All HRV metrics were computed exclusively from accepted RR intervals after artifact rejection and outlier correction.

As for EDA, each window was processed using NeuroKit2 (Makowski et al., 2021) to separate the slowly varying tonic component, or skin conductance level (SCL), from the faster phasic component containing skin conductance responses (SCRs). SCR events were automatically identified from the phasic component based on their onset, peak, and half-recovery points. The extracted features comprised mean EDA and mean SCL, representing the average overall and tonic skin conductance within each window, respectively; SCR peak count, defined as the number of detected phasic responses within a window; cumulative SCR amplitude (ΣSCR_*A*), calculated as the sum of the amplitudes of all detected responses; and mean SCR amplitude (SCR_*A*), calculated as their average amplitude. For each response, amplitude was measured as the increase in conductance from SCR onset to peak. Mean SCR rise time (SCRt_*rise*) represented the average interval from response onset to peak, whereas mean SCR recovery time (SCRt_*rec*) represented the average interval from the peak to the point at which the response decreased to half of its amplitude. These measures characterize the frequency, magnitude, and temporal dynamics of phasic sympathetic electrodermal activity (Posada-Quintero and Chon, 2020; Boucsein et al., 2012). It is worth noting that the introduction of phasic features was not intended to represent event-related SCRs linked to discrete stimulus onsets, but to characterize non-specific phasic activity throughout each scene. Similar studies have extracted features from continuous EDA signals without assigning responses to specific events (Veeranki et al., 2024; Kumar et al., 2023). Accordingly, these features were interpreted as window-level descriptors of the frequency, magnitude, and temporal morphology of sympathetic fluctuations.

Normalization was done through z-score, considering each participant's mean (x̄R) and standard deviation (s_*R*) at their 2-min baseline period. Thus, each non-normalized sample x was transformed into a normalized sample x′ according to Equation 2.

x′=(x-x̄R)sR(2)

Afterwards, to obtain an average value for each scene, the median of each brain region-frequency band combination and autonomic features was computed from normalized values, averaged within subjects. This median was used to create both the topoplots using Python (MNE, v1.5.1) and the multi-modal normalized heatmap using R (lattice, v4.5.1). Median instead of average across subjects was used to capture the sample's overall trend while being robust to outliers (Géron, 2019). Median Absolute Deviation (MAD; Equation 3) was calculated to measure the spread of the data, considering the absolute distance between each subject's data point (x_*i*) and the dataset's median (x~) (Géron, 2019).

MAD=median(|xi-x~|)(3)

2.6 Statistical analysis

Normalized values were used to detect differences between each fear scene and the baseline, as well as to determine whether normalized responses varied across the four fear scenes, using a two-sided exact permutation test and a Friedman test, respectively. p-values were Bonferroni-corrected to control for Type I error inflation due to multiple comparisons (Dunn, 1961).

Cohen's d is a measure of effect size that is standardized for the purpose of comparing the mean values of two groups (Lakens, 2013); effect size could be small (d = 0.2), medium (d = 0.5), or large (d = 0.8) (Cohen, 1988). Paired cohen's d (Rousselet et al., 2017) was calculated through Equation 4, where z¯ and s_*z* represent the mean and standard deviation of the differences between x₁ and x₂, which represent the average value for the fear scene and baseline, respectively.

dz=z¯sz;z=x1-x2(4)

Additionally, a more robust effect size was also calculated, considering the non-parametric paired Cliff's delta (Cliff, 1996), which estimates the probability that a random value from x₁ is greater than a random value from x₂, minus the reverse (Macbeth et al., 2011). Its ranges consist in negligible (0 <δ <0.14), small (0.14 <δ <0.33), medium (0.33 <δ <0.47), and large (δ > 0.47) (Romano et al., 2006). Cliff's δ was calculated as in Equation 5, where the cardinality symbol # indicates counting, x₁ the scores at the fear scene, x₂ the scores at the baseline, and N the number of pairs.

δ=#(x1>x2)-#(x1<x2)N(5)

Finally, a multi-modal hierarchical clustering was performed on overlapping subjects (n = 6; S2, S3, S4, S6, S8, S9) considering their subject-scene values, using Spearman distance and average linkage. Bootstrap resampling for cluster validation (pvclust, 1,000 iterations) was performed over these subject-scene values, and clusters were considered statistically supported at Approximately Unbiased (AU) p-values ≥95. Cophenetic correlation (cophenetic) and Baker's gamma (cor.dendlist) evaluated linkage effects across average, complete, ward.D2, and single methods. Finally, the Mantel statistic and quantile significance (mantel) tested cross-modal coupling using Spearman distance with 1,000 permutations. Spearman (Equation 6) was selected over Euclidean distance for its robustness to outliers (Boslaugh, 2012), where r_*s* of two features X and Y is calculated by their difference between ranks d, and the number of subjects n.

rs(X,Y)=1-6∑i=1ndi2n(n2-1)(6)

3 Results

For each fear scene, normalized neural (EEG; Table 1) and autonomic (HRV/EDA; Table 2), each table shows the median (MAD) value of each fear scene. Moreover, Tables 3, 4 show inferential statistics: Two-sided exact permutation test's p-value, paired Cohen's d and Cliff's δ.

Table 1. Median and median absolute deviation (MAD) of baseline-normalized EEG features.

The highest median value for each feature is shown in bold.

Table 2. Median and median absolute deviation (MAD) of baseline-normalized HRV and EDA features (winsorized at ± 3σ).

The highest median value for each feature is shown in bold.

Table 3. Inferential statistics on baseline-normalized EEG responses for each brain region and frequency band when compared to baseline.

Each fear scene was compared with baseline. Two-sided exact permutation test's Bonferroni-corrected p-value, paired Cohen's d and Cliff's δ. p <0.05, |d|>0.8 and |δ|>0.47 are in bold.

Table 4. Inferential statistics on baseline-normalized HRV and EDA responses when compared to baseline.

Each fear scene was compared with baseline.

Two-sided exact permutation test's Bonferroni-corrected p-value, paired Cohen's d and Cliff's δ. p <0.05, |d|>0.8 and |δ|>0.47 are in bold.

The normalized EEG topoplots are shown at Figure 2, and their inferential statistics at Figure 4, while the multi-modal z-score normalized heatmap with hierarchical clustering is shown at Figure 3.

3.1 Scene-based neurophysiological responses

Considering normalized, z-score EEG responses at Table 1 and inferential statistics at Table 3, there are consistent yet varied responses according to the scene in which the participant was situated. Considering the frontal region, on θ, there were no statistically significant z-score values, as it was slightly higher during the Conceptual scene (0.34, p = 1.16, d = 0.36, δ = 0.25), higher at Embodied (2.09, p = 0.06, d = 0.69, δ = 0.50) and Integrative (2.42, p = 0.10, d = 0.41, δ = 0.50), and highest at Discursive (4.12, p = 0.21, d = 0.54, δ = 0.50); α also do not showed statistically significant values, as it was slightly higher for all fear scenes (1.25, 0.18, 1.80, 0.89; p = 0.07, p = 1.30, p = 0.19, p = 0.07; d = 0.71, d = 0.36, d = 0.62, d = 0.50; δ = 0.75, δ = 0, δ = 0.50, δ = 0.75), with Discursive being the highest (1.80, p = 0.19, d = 0.62, δ = 0.50); and β had some statistically significant differences at Embodied, Discursive, and Integrative, as it was slightly higher during the Conceptual scene (0.08, p = 0.78, d = 0.37, δ = 0.00), higher at Integrative (2.80, p = 0.005, d = 0.54, δ = 0.75) and Embodied (3.19, p = 0.0003, d = 0.69, δ = 1.00), and highest at Discursive (11.87, p = 0.02, d = 1.02, δ = 0.75). Moreover, though not statistically significant across frequency bands, z-scored θ_Central power was slightly higher during the Integrative scene (0.81, p = 0.28, d = 0.65, δ = 0.50), higher at Discursive (3.07, p = 0.25, d = 0.55, δ = 0.50), and Embodied (3.25, p = 0.25, d = 0.55, δ = 0.25), and highest at Conceptual (5.23, p = 0.27, d = 0.87, δ = 0.25); similarly, α_Central was slightly higher during the Integrative scene (0.94, p = 0.38, d = 0.63, δ = 0.25), higher at Embodied (2.44, p = 0.30, d = 0.55, δ = 0.25), and Discursive (2.71, p = 0.28, d = 0.56, δ = 0.25), and highest at Conceptual (3.44, p = 0.31, d = 0.79, δ = 0.00); likewise, β_Central was slightly higher during the Integrative scene (0.23, p = 0.64, d = 0.60, δ = 0.00), higher at Discursive (2.53, p = 0.30, d = 0.42, δ = 0.25), and Embodied (3.28, p = 0.18, d = 0.65, δ = 0.50), and highest at Conceptual (6.15, p = 0.23, d = 0.65, δ = 0.50).

Among all EEG-derived indices, Fatigue was the most consistent and statistically significant across scenes, having its highest value at Discursive (1.79, p = 0.0003, d = 0.84, δ = 1.00), following Conceptual (1.77, p = 0.07, d = 0.43, δ = 0.75), Embodied (1.40, p = 0.0003, d = 0.60, δ = 1.00), and Integrative (1.21, p = 0.0003, d = 0.68, δ = 1.00); Excitement index also showed a similar pattern, with the highest value being at Embodied (5.94, p = 0.02, d = 0.54, δ = 0.75), followed by Integrative (3.60, p = 0.005, d = 0.53, δ = 0.75), Discursive (2.68, p = 0.0003, d = 0.49, δ = 1.00), and Conceptual (2.56, p = 0.21, d = 0.38, δ = 0.50); lastly, Engagement index had mixed responses, as both Embodied and Integrative had higher values (2.38, p = 0.04, d = 0.76, δ = 0.75; 2.12, p = 0.10, d = 0.69, δ = 0.50), while Discursive and Conceptual had lower values (0.94, p = 0.005, d = 0.53, δ = 0.75; 0.65, p = 0.35, d = 0.47, δ = 0.25).

Lastly, parietal and occipital regions exhibited mixed responses, with non-statistically significant values at θ_Parietal (–0.12, 0.04, 1.89, 0.83) and θ_Occipital (–0.21, –0.46, 3.81, 0.05), as well as decreasing z-score values at α on all scenes except for Discursive, for instance α_Parietal (–0.38, –0.39, 0.40, –0.20) and α_Occipital (–0.43, –0.47, 1.62, –0.40), with α_Parietal at Conceptual being the only statistically significant (–0.39, p = 0.02, d = −1.31, δ = −0.75); while β had increasing z-scores on all scenes, for instance β_Parietal (0.03, 0.54, 1.75, 0.56) and β_Occipital (0.77, 0.05, 3.30, 0.07), with β_Occipital at Discursive being the only statistically significant (3.30, p = 0.03, d = 0.63, δ = 0.50.

Moreover, channel-level normalized EEG responses are shown with topoplots (Figure 2) of θ, α, and β power for each fear scene (Embodied = E, Conceptual = C, Discursive = D, Integrative = I), as well as statistically significant differences to baseline with topoplots of Bonferroni-corrected exact permutation test's p-values (Figure 4). On θ power there were no statistically significant z-score values, as the Embodied scene had the highest activation at C4 (E = 2.01, C = 0.24, D = 1.60, I = 0.34) and FP2 (E = 2.32, C = 0.25, D = 0.35, I = 0.59); while the Discursive scene had the highest activation at O1 (E = –0.27, C = –0.33, D = 1.94, I = –0.21) and O2 (E = –0.07, C = –0.53, D = 1.60, I = 0.17), while the Integrative scene had the highest activation at P8 (E = –0.09, C = 0.13, D = 0.37, I = 1.33). Moreover, at the α power, there is a shared neural pattern of desynchronization at the parietal region, particularly, a significant decrease in P7 at Conceptual (E = –0.47, C = –0.43; p = 0.022, D = –0.37, I = –0.13) and a descriptively decrease in P8 (E = –0.45, C = –0.31, D = –0.04, I = –0.09); in addition to a synchronization in FP1 at Integrative (E = 0.53, C = –0.03, D = 1.00, I = 0.49; p = 0.027) and a descriptively increase in FP2 (E = 1.16, C = –0.03, D = 0.59, I = 0.49) channels, with Embodied and Discursive have non-statistically significant higher values at the C4 channel (E = 1.06, C = –0.09, D = 1.19, I = 0.07) but lower at C3 (E = –0.31, C = –0.06, D = –0.13, I = 0.18). Finally, regarding the β power, there is a significant synchronization pattern at Embodied, Discursive, and Integrative in FP1 (E = 3.13; p = 0.005, C = –0.05, D = 7.11; p = 0.022, I = –0.32; p = 0.005) and FP2 (E = 1.16; p = 0.001, C = 0.03, D = 6.87; p = 0.027, I = 2.64; p = 0.005) channels; in addition to significant increase in P8 power at Conceptual (E = –0.05, C = 1.62; p = 0.040, D = 1.73, I = 1.38) and O2 power at Discursive (E = 0.38, C = 0.37, D = 2.58; p = 0.027, I = –0.11).

Regarding heart rate and time-domain HRV metrics, none of the comparisons between the fear scenes and baseline remained statistically significant after Bonferroni correction (Table 4). Nevertheless, the descriptive statistics showed that the Discursive scene had the highest median values for SDNN, RMSSD, and pNN50 (Table 2). The largest differences relative to baseline, in terms of median, were observed for RMSSD (1.47, p = 0.31, d = 0.94, δ = 0.75) and pNN50 (2.60, p = 0.12, d = 1.20, δ = 0.75), both of which exhibited large effect sizes. BPM was descriptively higher than baseline during the Embodied, Conceptual, and Discursive scenes, with the largest standardized difference observed for the Embodied scene (2.25, p = 0.94, d = 0.62, δ = 0.50), although non-statistically significant. These findings should therefore be interpreted as descriptive patterns rather than statistically supported differences.

Regarding frequency-domain HRV metrics, only one comparison remained statistically significant after correction. Total power was lower during the Embodied scene (–0.66) than at baseline (p = 0.005, d = −1.10, δ = −0.75). Among the experimental scenes, Discursive showed the highest median values for HF power (0.90) and total power (0.50); however, neither difference relative to baseline was statistically significant (HF: p = 0.25, d = 0.98, δ = 0.50; total power: p = 1.000, d = 0.26, δ = 0.25). The LF/HF ratio was descriptively lower during all experimental scenes than at baseline, although none of these comparisons reached statistical significance. Thus, the frequency-domain results provided statistical support only for lower total power during Embodied, whereas the remaining differences represent descriptive patterns.

Regarding EDA, the tonic measures showed their highest median values during the Conceptual scene. Specifically, EDA appears to have increased from a median of 0 at baseline to 0.52 during Conceptual, although the comparison did not reach statistical significance after correction (p = 0.52, d = 0.82, δ = 0.50). SCL showed a similar pattern, increasing from a median of 0 at baseline to 0.53 during Conceptual (p = 0.52, d = 0.84, δ = 0.50). Thus, despite the large effect sizes, these differences should be interpreted as descriptive.

For the phasic EDA metrics, SCR_*P* was non-significantly lower during all experimental scenes than at baseline. The largest difference between medians was observed during Conceptual, decreasing from 0 at baseline to –0.62 (p = 1.00, d = −0.15, δ = −0.25). In contrast, both ∑SCR_*A* and SCR_*A* showed higher median values during all experimental scenes and reached their maxima during Discursive. The median of ∑SCR_*A* increased from 0 at baseline to 3.00 during Discursive, representing a statistically significant difference with large effect sizes (p = 0.005, d = 2.17, δ = 0.75). Similarly, the median SCR_*A* increased from 0 to 3.00 during Discursive (p = 0.005, d = 2.34, δ = 0.75). The largest differences between medians for the temporal SCR metrics were also observed during Discursive: SCR_{t_rise} increased from 0 at baseline to 1.70 (p = 0.56, d = 0.77, δ = 0.50), whereas SCR_{t_rec} increased from 0 to 0.52 (p = 1.00, d = 0.50, δ = 0.25). However, neither temporal comparison was statistically significant after correction. Overall, statistical support for the phasic EDA responses was limited to the higher ∑SCR_*A* and SCR_*A* values observed during the Discursive scene.

3.2 EEG-autonomic relationships

Further, a multi-modal normalized heatmap with a hierarchical cluster analysis is shown at Figure 3. Average linkage was selected for its highest cophenetic correlation (c = 0.70; complete c = 0.60; single c = 0.58; Ward (c = 0.53). Baker's γ indicated strong agreement within linkage families (average-single γ = 0.76; complete-Ward γ = 0.85), but weak agreement across families (γ= 0.19–0.25), indicating that structure is sensitive to linkage method. Additionally, a Mantel test on the Spearman correlation-based EEG and HRV distance matrices found a weak but significant cross-modal association (r = 0.082, p = 0.045).

Although several statistically significant clusters emerged after bootstrapping (AU≥95), they remained isolated across devices, having no significant higher-level connections (h>0.72, AU <95). Hierarchical clustering yielded five main branches: (1) a significant EDA cluster (EDA and SCL; h = 0.02, AU = 100), (2) a significant cross-modal HRV-EDA cluster (LF/HF and SCR_*P*; h = 0.45, AU = 99), (3) an EEG cluster, (4) a SCR cluster, and (5) a HRV cluster.

Within the EEG branch, frontal θ and α formed a significant cluster (h = 0.03, AU = 99) before merging significantly with β (h = 0.13, AU = 99). This group later merged significantly (h = 0.41, AU = 96) with a non-significant EEG index cluster, composed of a non-significant Engagement-Excitement group (h = 0.14, AU = 92) and Fatigue (h = 0.38, AU = 83). At h = 0.72, this entire frontal and EEG index cluster merged significantly (AU = 99) with a compound parietal-occipital-central cluster. Here, a non-significant parietal sub-cluster (θ and α at h = 0.15, AU = 83, joined non-significantly by β at h = 0.28, AU = 92) joined a significant occipital-central sub-cluster at h = 0.45 (AU = 99). That sub-cluster was non-significantly formed at h = 0.35 (AU = 93) by combining a non-significant occipital group (α and β at h = 0.02, AU = 88, joined significantly by θ at h = 0.04, AU = 100) with a non-significant central group (θ and α at h = 0.03, AU = 93, joined significantly by β at h = 0.09, AU = 100).

For SCR features, a non-significant sub-cluster of ΣSCR_*A* and SCR_*A* (h = 0.30, AU = 94) joined a significant SCR_trise-SCR_trec cluster (h = 0.34, AU = 95) to form a significant SCR branch at h = 0.60 (AU = 95). Finally, within HRV features, a non-significant HF-Total P group (h = 0.12, AU = 91) joined LF significantly (h = 0.24, AU = 95). This branch merged non-significantly at h = 0.50 (AU = 93) with a secondary branch formed by a non-significant RMSSD-pNN50 group (h = 0.22, AU = 82) and SDNN (h = 0.34, AU = 90). Both sub-clusters ultimately merged significantly with BPM at h = 0.63 (AU = 99).

4 Discussion

Participants exhibited shared neural and autonomic patterns across all fear scenes compared to baseline, identifying potential markers of fear. Additionally, unique responses at each scene also showed contextual modulation. Finally, multi-modal integration was essential for understanding the psychophysiological interplay between neural and autonomic dynamics during these fear scenes.

4.1 Neural and autonomic signatures of fear

The most consistent scene-invariant trend was the statistically significant increased activation of β_Frontal for all scenes except for the Conceptual scene (Tables 1, 3), also seen at the channel-level considering FP1 and FP2 (Figures 2, 4). Frontal regions are critical on fear-related neural circuity (Åhs et al., 2015; Picó-Pérez et al., 2019; Milad and Quirk, 2012; Rezaei et al., 2025; Kim et al., 2019; DeLaRosa et al., 2014; Apicella et al., 2023; Demaree et al., 2005), as they play a central role in fear and cognitive control (Kim et al., 2019; Demaree et al., 2005), and β in particular has been associated with increased attentional demands (Kim et al., 2019), perceptual processing (Okazaki et al., 2008), action execution (Lamarche and Ogilvie, 1997), and fear-related cognition (Apicella et al., 2023). Moreover, there was a descriptive increase of θ_Central and θ_Frontal across scenes, which has been related with increased decision-making and planning (DeLaRosa et al., 2014), as well as mental effort (Candela-Leal et al., 2025b; Villafaina et al., 2021b). Furthermore, α_Parietal showed a significant decrease at the Discursive scene, as well as descriptive decrease at Embodied and Integrative scenes, this region has been associated with emotional processing (Ugalde-Canitrot et al., 2025; Alemayehu et al., 1995; Yaşgüçlükal et al., 2021) and visual spatial attention (Rihs et al., 2007; Sauseng et al., 2005; Mann et al., 1996; Villafaina et al., 2021a), with its α downregulation linked to immersive virtual environments (Botrel et al., 2025; Magosso et al., 2019). Consequently, given Fatigue (θ/α), Excitement (β/α), and Engagement (β/(θ+α)) indices definition, and that both θ_Central and β_Frontal increased, as well as that α_Parietal decreased, EEG-derived indices significantly increased across all scenes except for the Conceptual scene, and Engagement for the Integrative scene. This was expected as the increased visual and auditory stimulation often increases mental weariness and alertness (Ramírez-Moreno et al., 2021; Pfurtscheller and Lopes da Silva, 1999; Pope et al., 1995; Candela et al., 2026; Romo-De León et al., 2026). Ultimately, the statistically significant modulations of α and β bands in are in line with earlier fear-related processing findings that considered θ and β bands (Rezaei et al., 2025), as well as α, β, and γ bands (Russo et al., 2024).

Regarding autonomic responses (Tables 2, 4), BPM was descriptively higher than baseline during the Embodied, Conceptual, and Discursive scenes, although none of these differences remained statistically significant after correction. This pattern may be compatible with increased autonomic arousal, effort, or threat-related activation during immersive tasks (Stiedl and Hager, 2017; Luong and Holz, 2022; Malińska et al., 2015). This is to mean that the observed behavior cannot be attributed exclusively to fear. The experimental scenes differed from the resting baseline in motor activity, cognitive engagement, verbal demands, and audiovisual stimulation. These factors may independently affect BPM and may also introduce variation or artifacts into wrist-recorded signals (Nelson et al., 2020; Menghini et al., 2019; Aisy Al et al., 2024; Shaffer and Ginsberg, 2017; Luong and Holz, 2022; Malińska et al., 2015; Lebamovski et al., 2025). Accordingly, the observed BPM pattern could be more appropriately interpreted as reflecting the combined demands of each scene rather than a fear-specific fight-or-flight response.

For EDA, the significantly higher ∑SCR_*A* and SCR_*A* values during the Discursive scene suggest increased phasic sympathetic activity under this condition. This interpretation is consistent with studies reporting larger or more frequent SCRs during fear, stress, and immersive experiences (Ghiasi et al., 2020; Yu et al., 2024; Kreibig et al., 2007; Luong and Holz, 2022). Nevertheless, these window-level metrics represent non-specific phasic activity and cannot be linked to individual fear-related events. Moreover, the remaining EDA features, including the relatively lower SCR_*P* values and the higher SCR t_rise and SCR t_rec values during some scenes, did not remain statistically significant after correction and should therefore be regarded as exploratory findings. An interesting pattern is that of the relatively higher pNN50 and RMSSD values during the Discursive scene and the generally lower LF/HF ratios, as they may hint parasympathetic predominance (which contradicts the observed phasic behavior) (Thomas et al., 2019; Shaffer and Ginsberg, 2017; Kleiger et al., 2005), however, this should not be interpreted as conclusive evidence given the lack of statistical significance. Moreover, these indices are sensitive to respiratory patterns, which were not monitored or controlled during the experiment (Schipke et al., 1999; Grossman and Taylor, 2007). Verbal activity, movement, task engagement, and the characteristics of the semi-immersive environment may therefore have contributed to the apparently mixed autonomic patterns (Nelson et al., 2020; Menghini et al., 2019; Aisy Al et al., 2024; Shaffer and Ginsberg, 2017; Luong and Holz, 2022; Malińska et al., 2015; Lebamovski et al., 2025). Overall, the findings suggest a complex autonomic response to the scenes, but the present design does not allow the effects of fear to be separated from those of respiration, behavioral demands, sensory stimulation, and measurement-related variation. Further studies incorporating respiratory and movement monitoring, together with appropriately matched active control conditions, are needed to clarify these contributions.

4.2 Contextual modulation of fear

Even though there were shared autonomic responses and neural activity across scenes, some had distinct modulations; thus, characterizing participants' responses to each scenario. While permutation tests revealed statistical significant differences between each scene and the baseline, a Friedman test testing for overall differences among the four scenes themselves was non-significant across EEG, HRV, and EDA (p>0.05). Inter-scene differences are thus treated as descriptive rather than confirmatory.

Even though the Embodied, Discursive, and Integrative scenes showed statistically significant increases in β_Frontal, β_FP1, and β_FP2, the Conceptual scene did not. Interestingly, the Conceptual was the only scene to exhibit an statistically significant decrease of α_Parietal, and the channel-level topoplot (Figure 4) revealed a left-lateralized decrease at α_P7, aligning with previously reported emotional lateralization effects in immersive fear paradigms, which reflect diverse cognitive demands during fear processing (DeLaRosa et al., 2014; El Basbasse et al., 2023). For instance, β_Frontal is related to perception (Kim et al., 2019; Okazaki et al., 2008) and action execution of fear (Lamarche and Ogilvie, 1997; Apicella et al., 2023), while α_Parietal to visual spatial attention (Rihs et al., 2007; Sauseng et al., 2005; Mann et al., 1996; Villafaina et al., 2021a) and threat salience (Chien et al., 2017; Codispoti et al., 2023; Arana et al., 2022). Moreover, despite similar Fatigue and Excitement indices across scenes, there was reduced Engagement for Conceptual and Discursive scenes when compared to other scenes, and given indices definition, these differences are derived from increased θ_Frontal power, related with enhanced cognitive planning processes (DeLaRosa et al., 2014). Another interesting finding was the descriptive increase of β_Parietal and β_P8 activity during Discursive and Integrative fear, as well as the statistically significant increase of β_P8 during Conceptual fear, which has been associated with conceptual-level fear generalization (and Morey et al., 2020), belief updating (Hein and Herrojo Ruiz, 2022), and memory-level threat processing (Wannemueller et al., 2024); while decreases have been linked to attention to threat and motor preparation (Sperl et al., 2021; Pirazzini et al., 2023); and right-hemispheric predominance has been related to negative emotion processing (El Basbasse et al., 2023; Tomarken et al., 1990; Prete et al., 2015; Ahern and Schwartz, 1979). These results suggest distinct cognitive demands across conditions: Embodied elicited a bottom-up sensory processing and action execution, Conceptual engaged top-down belief updating and threat salience, while Discursive and Integrative scenes involved a blend of both.

On the other hand, the autonomic pattern observed during Discursive fear may suggest context-dependent modulation, as this scene showed some of the largest descriptive differences in HRV and EDA metrics relative to baseline. Although the higher BPM and EDA responses observed during some scenes are compatible with sympathetic activation associated with fear exposure and threat-related arousal (Stiedl and Hager, 2017; Luong and Holz, 2022; Malińska et al., 2015), they cannot be attributed exclusively to fear. Movement, cognitive engagement, verbal demands, and other characteristics of the experimental tasks may independently influence autonomic activity (Nelson et al., 2020; Menghini et al., 2019; Aisy Al et al., 2024). During Discursive fear, ∑SCR_*A* and SCR_*A* were significantly higher than at baseline, suggesting increased phasic sympathetic activity under this condition. However, the Discursive scene also showed descriptively higher RMSSD, pNN50, and HF power, although these differences were not statistically supported. Considering the nature of the task, this apparently mixed autonomic pattern may partly reflect the verbal and respiratory demands of emphatic speech, as speech production can simultaneously increase sympathetic activation and modify vagally mediated HRV through respiration-related changes during vocalization (Dodo and Hashimoto, 2019). Because HF power, RMSSD, and pNN50 are strongly influenced by respiration (Schipke et al., 1999; Grossman and Taylor, 2007), these descriptive differences may reflect speech-related respiratory modulation rather than parasympathetic activation attributable to fear. Moreover, verbalization itself has been associated with regulatory effects during unpleasant emotional processing, including reduced sympathetic dominance and increased parasympathetic-related activity (Rantanen et al., 2021). Thus, the autonomic pattern observed during Discursive fear may reflect cognitively mediated and context-dependent processing together with the specific behavioral demands of this scene, rather than a purely reflexive threat response. However, the Friedman tests identified no significant differences among the four fear scenes. Therefore, the present results do not provide strong enough statistical support for scene-specific autonomic modulation, and the comparatively pronounced pattern observed during Discursive fear should be regarded as descriptive and hypothesis-generating.

4.3 Multi-modal integration

Multi-modal hierarchical clustering with bootstrap validation on overlapping subjects revealed a modality-specific organization rather than broad brain-autonomic coupling. Five main branches aligned with signal type: EDA (AU = 100), SCR (AU = 95), HRV (AU = 99), and EEG, with neural and autonomic branches remaining statistically independent (h>0.72, AU <95). The only validated cross-modal link (HRV-EDA) joined LF/HF with SCR_*P* (h = 0.45, AU = 99), consistent with shared parasympathetic/sympathetic influence on autonomic arousal (Shaffer and Ginsberg, 2017; Thomas et al., 2019). Within EEG, the frontal region (θ, α, β) formed a stable cluster (AU = 99), consistent with attentional engagement during sustained threat processing (Rihs et al., 2007; Sauseng et al., 2005); however, deeper merges with parietal-occipital-central structures were non-significant (AU <95). Finally, a weak yet significant global EEG-autonomic Mantel association (r = 0.082, p = 0.045) indicated shared variance without the cross-modal clustering required to confirm strong brain-body coupling.

Concurrent brain and autonomic responses highlight the value of multi-modal approaches for enabling a more comprehensive characterization of fear-related processing across brain and body domains, providing greater robustness than unimodal analyses (Ramírez-Moreno et al., 2026; Cos et al., 2023) and supporting its potential utility for interactive and immersive learning environments (Bălan et al., 2020a,b; Freeman et al., 2018).

4.4 Limitations

Although consistent neural and autonomic responses were observed across fear-related scenes, together with context-dependent modulation of fear neurophysiology measured through EEG and HRV/EDA, some limitations should be acknowledged.

First, the relatively small sample size (N = 10), along with the reduced effective sample sizes for inter-scene comparisons (N = 8) and multi-modal hierarchical clustering (N = 6), limits statistical generalization; however, multi-modal neurophysiological studies that integrate EEG and peripheral signals are commonly constrained to small cohorts due to device synchronization, pre-processing, and experimental complexity (Ramírez-Moreno et al., 2026), reporting sample sizes of 1 (Pineda-Hernández, 2022; Fuentes-García et al., 2022), 10 (Botrel et al., 2025), 11 (Villafaina et al., 2021b), 12 (Romo-De León et al., 2024), 13 (Olivas et al., 2021), and 14 (Ayuso-Moreno et al., 2020; Villafaina et al., 2019; Yu et al., 2024). Consequently, the present exploratory study emphasized effect sizes and overall response patterns, supported by Cohen's d and Cliff's delta, rather than population-level inference. Future larger cohorts studies could validate both shared and scene-specific findings.

While scenes' variations in task demands, movement, and audiovisual stimulation were inherent to the design's ecological validity, they introduce potential confounding factors. Although responses were normalized to the neutral baseline, the absence of non-fear controls means we could not isolate fear from task or environment-driven effects, such as α band downregulation and increased engagement associated with immersive environments (Botrel et al., 2025; Candela et al., 2026), or cumulative carry-over effects from the fixed order constructive design. Nevertheless, ~1-min washout periods were used to explain the next scene, interrupt immersion, and mitigate physiological carry-over, future studies should also employ continuous recordings to quantify if these intervals returned signals to baseline. In addition, respiration was not controlled and may have influenced autonomic measures such as HF power, RMSSD, and pNN50 (Schipke et al., 1999; Grossman and Taylor, 2007). Moreover, given that the audiovisual environment was adaptively modulated using real-time EEG-derived VAD predictions (Ortiz-Espinoza et al., 2025; Blanco-Ríos et al., 2024), EEG served as a control signal and outcome measure. Hence, the observed EEG differences should be interpreted within the closed-loop paradigm rather than only stimulus-driven effects. Accordingly, while not pilot-testing Scene 1 stimuli for fear induction limits inter-studies direct comparisons, it aligns with the paradigm's ecological focus on adaptive, emotion-modulated interactions.

Besides this, the autonomic findings should not be interpreted as responses exclusively attributable to fear. The resting baseline and experimental scenes were not matched for motor activity, cognitive engagement or audiovisual stimulation, and these factors were not experimentally controlled or managed during data analysis. Movement, in particular, may introduce artifacts into wrist-recorded data, further complicating the distinction between physiological activation and measurement-related variation (Shaffer and Ginsberg, 2017; Boucsein et al., 2012). Consequently, comparisons with baseline should be understood as reflecting the combined autonomic response to each scene rather than a fear-specific effect. Future studies may incorporate respiratory and movement monitoring and include active control conditions matched for the principal motor, cognitive, verbal and sensory characteristics of each scene. Additionally, they may incorporate motion artifact rejection using accelerometry data as a reference signal (e.g., adaptive filtering).

5 Conclusions

The present exploratory study identified neural and autonomic signatures of fear during concurrent EEG and peripheral recordings in a semi-immersive environment across multidimensional fear-related scenarios, including Embodied, Conceptual, Discursive, and Integrative.

Shared fear-related responses included frontal and parietal brain activity, increased EEG-derived indices, and elevated BPM and SCR_*A* measures, while context-dependent modulation was reflected by scene-specific brain lateralization patterns, with increased β_Frontal activity during Embodied fear and decreased α_Parietal activity during Conceptual fear, consistent with bottom-up and top-down processing, respectively. Similarly, Discursive fear showed mixed autonomic responses with speech-related respiratory modulation and cognitively mediated rather than reflexive threat processing.

Moreover, multi-modal clustering revealed that fear responses are mainly organized within, rather than across, neural and autonomic systems: features formed robust modality-specific clusters, and while a global EEG-autonomic Mantel association was weak but significant, the only validated cross-modal cluster linked HRV (LF/HF) and EDA (SCR_*P*). Despite sample size limitations, these findings highlight the value of multi-modal immersive paradigms for capturing distinct and shared dimensions of fear physiology.

Generative AI statement

The author(s) declared that Generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Atsushi Noritake, University of Toyama, Japan

Reviewed by: Luciano Cavallino, National Scientific and Technical Research Council (CONICET), Argentina

Shaoting Zeng, Beijing University of Technology, China