Assessing Autonomic Function from Electrodermal Activity and Heart Rate Variability During Cold-Pressor Test and Emotional Challenge

Standard functional assessment of autonomic nervous system (ANS) activity on cardiovascular control relies on spectral analysis of heart rate variability (HRV) series. However, difficulties in obtaining a reliable measure of sympathetic activity from HRV spectra limits the exploitation of sympatho-vagal metrics. On the other hand, measures of electrodermal activity (EDA) have been demonstrated to provide a reliable quantifier of sympathetic dynamics. In this study we propose novel indices of phasic autonomic regulation mechanisms by combining HRV and EDA correlates and thoroughly investigating their time-varying dynamics. HRV and EDA series were gathered from 26 healthy subjects during a cold-pressor test and emotional stimuli. Instantaneous linear and nonlinear (bispectral) estimates of vagal dynamics were obtained from HRV through inhomogeneous point-process models, and combined with a sensitive maker of sympathetic tone from EDA spectral power. A wavelet decomposition analysis was applied to estimate phasic components of the proposed sympatho-vagal indices. Results show significant statistical differences for the proposed indices between the cold-pressor elicitation and previous resting state. Furthermore, an accuracy of 73.08% was achieved for the automatic emotional valence recognition. The proposed nonlinear processing of phasic ANS markers brings novel insights on autonomic functioning that can be exploited in the field of affective computing and psychophysiology.

Assessing Autonomic function from electrodermal Activity and Heart Rate Variability During cold-pressor test and emotional challenge Shadi Ghiasi 1* , Alberto Greco 1 , Riccardo Barbieri 2 , enzo pasquale Scilingo 1 & Gaetano Valenza 1 Standard functional assessment of autonomic nervous system (AnS) activity on cardiovascular control relies on spectral analysis of heart rate variability (HRV) series. However, difficulties in obtaining a reliable measure of sympathetic activity from HRV spectra limits the exploitation of sympatho-vagal metrics. on the other hand, measures of electrodermal activity (eDA) have been demonstrated to provide a reliable quantifier of sympathetic dynamics. In this study we propose novel indices of phasic autonomic regulation mechanisms by combining HRV and eDA correlates and thoroughly investigating their time-varying dynamics. HRV and EDA series were gathered from 26 healthy subjects during a cold-pressor test and emotional stimuli. instantaneous linear and nonlinear (bispectral) estimates of vagal dynamics were obtained from HRV through inhomogeneous point-process models, and combined with a sensitive maker of sympathetic tone from eDA spectral power. A wavelet decomposition analysis was applied to estimate phasic components of the proposed sympatho-vagal indices. Results show significant statistical differences for the proposed indices between the cold-pressor elicitation and previous resting state. Furthermore, an accuracy of 73.08% was achieved for the automatic emotional valence recognition. the proposed nonlinear processing of phasic AnS markers brings novel insights on autonomic functioning that can be exploited in the field of affective computing and psychophysiology. The sympathetic nervous system (SNS) plays a fundamental role in driving cardiovascular control dynamics together with the other branch of the Autonomic Nervous System (ANS), the parasympathetic nervous system 1,2 , and includes the regulation of blood pressure, respiration, and heart rate 1,3,4 . Many sympathetic neurons are part of the central nervous system (CNS) having origins in the thoracic and lumbar regions (T1-L2) in the spinal cord. Most of the sympathetic activity control on cardiovascular functioning is performed through noradrenergic neurons. In case of a drop in arterial pressure during situations such as anger, fear or excitement, the arterial baroreflex senses these changes in blood pressure through the baroreceptors, and the ANS acts to increase the vasoconstriction. As a result, the afferent input to the central autonomic nuclei is decreased, which in turn results in an increase in the sympathetic neural outflow. The resulting sympathetic influences in conjunction with the parasympathetic ones leads to an increase in heart rate, oxygen intake, and the blood supply to the heart and skeletal muscles. The release of epinephrine (adrenaline) and norepinephrine (noradrenaline), as well as the stimulation of sweat glands are also the results of a sympathetic nervous system elicitation 2,5 . The sympathetic vasodilator nerves in the human skin results from a cholinergic co-transmission, which releases acetylcholine as a primary neurotransmitter. The acetylcholine binds to muscarinic receptors on sweat glands to activate sweating, defining a sympathetic skin response 6 . An early quantification of the sympathetic activity has been performed through measurement of plasma or urinary norepinephrine. However, these indices may not be reliable due to their low sensitivity and non specific localization of the sympathetic response 1 . While a standard sympathetic

Materials and Methods
Subjects recruitment, experimental protocol, and acquisition set-up. Twenty-six right-handed young healthy volunteers (8 females) from the University of Pisa participated in the study. Subjects did not have any history of neurological and cardiovascular diseases, or alcoholic or smoking habits. They were asked to avoid coffee and alcohol, as well as strenuous exercise, at least 2 hours before the laboratory visit [39][40][41] . Subjects aged between 21-32 (26 ± 3) with 60-100 kilograms of weight and 160-175 centimeters of height. All participants were screened using Patient Health Questionnaire (PHQ) 42 . This test intends to exclude the subjects with mental disorders for the experiment since the population under study is healthy control. Therefore, only participants with a PHQ score lower than 5 were included in the study. Valsalva maneuver and breathing holding was forbidden throughout the study. To ensure the hemodynamic stabilization, before the experiment participants were asked to sit in a comfortable chair, while watching a black screen. The experiment comprised two main phases: Cold-Pressor Test (CPT) and affective elicitation phase. The protocol timeline was as follows (see also Fig. 1): • 4 minutes of resting state • up to 3 minutes of CPT • 3 minutes of resting state • 4 minutes and half of emotional stimuli through videos • 4 minutes of recovery.
The presentation order of CPT and emotional stimuli was randomized among subjects. Concerning the CPT, participants were asked to submerge their left hand (i.e, the non-dominant hand) up to wrist into a tank filled of ice and water (crushed ice) with the temperature of 0-4 degrees centigrade for a period of 3 minutes 43 . The choice of 3 minutes is consistent with the average pain threshold of healthy subjects 44 . In case a subject could not tolerate this kind of elicitation, they were quickly moved to the next resting phase.
The emotional elicitation phase foresaw video-clips using a video-projector. The well-known Circumplex Model of Affect (CMA) 45 was used for the emotional modeling through its valence and arousal dimensions. Valence refers to the pleasantness of the emotional perception, whereas arousal represents the strength of such a perception. In our study, the three categories of video-clips were associated with the following sessions: (1) high arousal and high positive valence, (2) neutral arousal/valence, (3) high arousal and high negative valence. The three categories have the same duration of 1 minute and 30 seconds and the presentation of the clips was randomized among subjects. The experiment ended by asking subjects for a complete self-evaluation of the stimuli through the visual analog scale questionnaire (VAS) 46 . Throughout the experiment, continuous electrocardiogram (ECG) and EDA signals were acquired using the BIOPAC MP35 device with a sampling rate of 500 Hz for further analyses 47 . The dedicated modules, ECG100C Electrocardiogram Amplifier and EDA100C Electrodermal Activity Amplifier from BIOPAC inc. were used to acquire the ECG signal and the electrodermal activity, respectively. The ECG100C ECG Amplifier has a noise voltage of 0.1 V (rms) in the bandwidth 0.05-35 Hz, a CMRR of 110 dB, a 16-bit ADC, and an amplitude resolution of 0.4 mV to record the D2 lead ECG signal. Pregelled Ag/ AgCl electrodes for the ECG acquisition were placed following the Einthoven triangle configuration. EDA sensors were placed on the distal phalanx of the second and third finger of the dominant hand, imposing a DC voltage of 0.5V.
The RR series, defined as the interval between two successive QRS complexes, were derived by applying the automatic QRS complex detection algorithm on acquired ECG signals 48 . The technical artifacts due to the errors in R-peak detection algorithm were removed by applying cubic spline interpolation method 49 . The implementation was performed using Matlab software (R2017b version) and the Kubios HRV software 50 . The experimental procedure was approved by the "Comitato Etico Regionale per la Sperimentazione Clinica della Regione Toscana", section "Area Vasta Nord Ovest" -Protocol n. 7803, Registry number 1072, approved on 18 Jan 2018. The recordings were carried out in agreement with the Declaration of Helsinki. Written informed consent was obtained from all subjects.
instantaneous estimates of vagal activity from inhomogeneous point-processes for heartbeat dynamics. Full details in this methodological framework are reported in 20 , referring to the recently proposed nonlinear parametric point-process model continuously characterizing the probabilistic properties of heartbeat events.
Briefly, assuming the history dependence, an inverse-Gaussian probability density function (f(t|H t , ξ(t))), parametrized in its first-order moment (i.e., the instantaneous mean μ RR (t, H t , ξ(t))) and the shape parameter , characterizes each heartbeat event as a function of the history H t = (u j , RR j , RR j−1 , . . . , RR j−M+1 ), with , which are estimated through a maximum log-likelihood procedure with Newton-Raphson algorithm 20 . Note also that this formulation is used to pre-process all heartbeat data to identify and eventually correct errors including peak mis-detection and ectopic beats 51 .
The first-order moment is defined as a nonlinear Volterra-Wiener autoregressive model: is the output of the Laguerre filters before time t and is the i th -order discrete time orthonormal Laguerre function, with (n ≥ 0) and α the discrete-time Laguerre parameter 20 . Kolmogorov-Smirnov (KS) test and associated KS statistics are employed to calculate the model goodness-of-fit 20 . The independence of the model-transfomed intervals are tested with the use of Auto-correlation plots 20 . Through this framework, three levels of quantitative tools defined in the time, spectral, and bispectral domains can be obtained. Time domain estimates refer to the first-order (μ RR (t, H t , ξ(t))) and second-order (σ RR (t, H t , ξ(t))) moments of the underlying probability function f(t|H t , ξ(t)) characterizing the heartbeat series at each moment in time 20 . Spectral quantifiers come from the linear power spectrum estimation Q(f, t, H t , ξ(t)), which depends on the linear terms g 1 (i, t) and related transformation from the Laguerre space to the Wiener-Volterra input-output terms 20 . Once estimated, the continuous Q(f, t, H t , ξ(t)) can be integrated within the low frequency (LF = 0.05-0.15 Hz) and high frequency (HF = 0.15-0.5 Hz) bands to derive instantaneous sympatho-vagal and vagal estimates, respectively. On the bispectral domain, the bispectrum |Bis(f 1 , f 2 , t)| is defined from the linear and nonlinear input-output Wiener-Volterra terms, and can be integrated in the LF and HF bands to derive the following quantifiers: = .
. . Note that the bispectral quantifier HH has been demonstrated to provide estimates of vagal autonomic outflow 52 . All the introduced estimates are calculated with a 5 ms time resolution. More details about the formulation of the bispectrum can be found in 53 . instantaneous estimates of sympathetic activity from eDA. The EDA signal reflects changes in the skin conductance induced by the sweat gland activity, which is directly controlled by the sympathetic branch of the ANS. Therefore, it is possible to rely on EDA signal processing to quantify the sympathetic dynamics through a proper functional analysis in the frequency domain 7,15 .
Specifically, a normalized EDA series, as a result of a Z-score transformation, is downsampled to 50 Hz. Afterwards, a time-frequency representation of the series is obtained through short-time Fourier transform and Welch periodogram, with a Blackman window of 60s and 59s overlap. Hence, a time-varying estimate of sympathetic outflow EDA Symp(k) in k discrete time samples, within the interval (0, T), is obtained by integrating the time-frequency plane within the (0.045-0.25 Hz) band 7,15 .
Previously, a convex optimization modeling approach has been proposed for the decomposition of the skin conductance measure of EDA 54 . As the result of this decomposition, low and high frequency components of the skin conductance response (SCR) known as tonic and phasic activities, respectively, are separated. Mainly the SCR has been quantified by calculating the number of significant (above threshold) phasic driver peaks (EDA.nSCR), the sum of SCR amplitudes with respect to significant peaks (EDA.SumAmpSCR), the maximum value of phasic activity (EDA.PhasicMax). The quantification of the tonic component is done through the mean tonic activity computation (EDA.Tonic) as well. www.nature.com/scientificreports www.nature.com/scientificreports/ proposed instantaneous indices of sympatho-vagal balance and related phasic dynamics. Aiming to derive reliable sympatho-vagal balance estimators, we combined the aforementioned EDA Symp (k), taken as a consistent measure of the sympathetic tone, and point-process derived measures from HRV linked to vagal activity, namely spectral power HF and bispectral powers LL, LH, and HH. Specifically, the following sympatho-vagal measures are defined: The derivation of sympatho-vagal measures is summarized in Fig. 2. We take a step further characterizing sympatho-vagal balance by (i) decomposing their time-varying dynamics into two components (1) the tonic activity representing a slow trend dynamic and (2) superimposed phasic component known to reflect high frequency oscillations of the time series and (ii) quantifying phasic dynamics using the median and area under the curve (AUC) metrics. Tonic and phasic decomposition of time-varying bispectral estimates LL, LH, and HH is further performed. An exemplary decomposition is shown in Fig. 3. Using a wavelet decomposition procedure 55 , tonic (S ton ) and phasic (S ph ) dynamics were derived from the first level of approximation coefficients and fifth level of detail coefficients, respectively, with a Daubechie5 function as the mother wavelet 55,56 . Classification for automatic valence recognition. Features were grouped into six categories. The first category F 1 is from HRV measures, including instantaneous point-process estimates defined in the time and frequency domains and standard HRV time domain measures including Root Mean Square of the Successive Differences (RMSSD) and normalized mean number of times an hour in which the change in successive normal sinus (NN) intervals exceeds 50 ms (pNN50).   Table 1.
Time-varying dynamics were averaged within a 90s time window, which corresponds to the duration of either the neutral, pleasant and unpleasant elicitation. In order to minimize the inter-subject variability, for each subject the feature set related to a negative or positive stimulus was normalized by dividing for the value computed during a neutral session. These features were used to automatically classify emotional states associated with negative and positive elicitations through a standard nonlinear Support Vector Machine (SVM) with recursive feature elimination (RFE), which was validated through a Leave-One-Subject-Out (LOSO) procedure. RFE is a wrapper feature selection method that chooses a subset of features with size r among m features (r < m) that optimizes the performance of the SVM classifier 57,58 . The method is based on a backward sequential selection. Particularly, the feature that has the least impact on the SVM weight-vector norm is eliminated for each repetition. Therefore, the features are ranked and the SVM classification is repeated m times while the last ranked features are removed at each iteration. In this study, a recently developed nonlinear SVM-RFE which is an improved version of the previous RFE algorithm is implemented 59 . This developed method tackles the issue of the correlation bias of the previous method which makes it reliable if the input feature set contains highly correlated features.
The LOSO validation scheme foresees that, at each iteration, the SVM-RFE model is trained using feature samples from N − 1 subjects (where N is the total number of subjects), and tested on the feature samples from the left-out subject. This is iterated N times. Final feature ranking was obtained by summation of the ranks from each iteration.
Classification results are reported in terms of confusion matrix associated with the feature subset giving highest accuracy. Note that diagonal values of this matrix represent sensitivity and specificity of classification, which are defined as the proportion of actual positives and negative valence elicitations that were correctly identified, respectively. The recognition was performed using MATLAB v8.4 and an additional toolbox for classification 60 . The processing chain of the classification for valence recognition is depicted in Fig. 4.

experimental Results
We first report on a comprehensive characterization in the statistical sense of the proposed features performance by comparing resting vs. CPT phases. Given the non-Gaussianity of samples, as revealed by p < 0.05 from Kolmogorov-Smirnov tests with null hypothesis of data Normality, p-values for the Rest vs. CPT comparison are from Wilcoxon non-parametric tests for paired data, with null hypothesis of equal medians between samples. These results aim at ensuring that the phasic dynamics of the proposed linear and nonlinear sympatho-vagal measures carry significant information on ANS changes. Afterwards, we exploit this knowledge in recognizing positive vs. negative valence during the emotional video elicitations by applying the aforementioned SVM-RFE algorithm.  Table 2, where the feature time-varying evolution is expressed as Median ± . Among the F 1 feature set (HRV time and frequency domain features), only μ RR showed significant differences between rest and CPT for all the selected time windows. Significant differences were found also in HF pp for the 30s and 180s windows. In the F 2 feature set (HRV bispectral features), significant differences were found for LL pp at all time windows, as well as for LH pp taking estimates from 180s CPT. Remarkably, significant differences were found for the proposed features in F 3 and F 5 , combining sympathetic EDA with linear and nonlinear vagal dynamics from HRV, but S LL in the 180s time window. Concerning features from phasic dynamics in F 6 , significant p-values were found for S LL ph for the 90s and 180s time windows. However, calculating the AUC of such phasic dynamics, bispectral features showed differences at 90s, whereas the ones from nonlinear sympatho-vagal balance were all different except for S HH ph estimated at 30s. www.nature.com/scientificreports www.nature.com/scientificreports/ Valence recognition results. The recognition accuracy as a function of the number of features, ranked through the SVM-RFE procedure, is shown in Fig. 6, whereas the confusion matrix associated with the best classification accuracy is reported in Table 3. This was achieved using four features, namely S LH ph , LF pp , LH ph , and S HH ph , with a balanced accuracy of 73.08% resulting in 75% and 71.43% of positive predictive value and negative predictive values, respectively.
Note that such most informative features come from spectral power in the LF, as well as from the proposed phasic dynamics of bispectral sympatho-vagal features. The complete feature set used for the classification, ranked through the SVM-RFE procedure, is reported in Table 4. As a comparison the classification results using other feature sets including only the indices extracted from HRV (F 1 and F 2 ) and the standard EDA indices (F 4 ) are also reported (Tables 5 and 6, respectively).

Discussion and conclusion
In this study we investigated novel ANS metrics linked to sympatho-vagal dynamics by combining EDA and HRV spectral powers, as well as combining EDA spectral power with HRV bispectral power. To this end, while sympathetic estimates were gathered from the electrodermal activity (EDA), vagal estimates were from HRV series. Features from phasic activity series of HRV bispectral power and combined EDA spectra were investigated as well. Preliminary results of this study can be found in 22,23 .
More in detail, we combined features from EDA power spectra, namely EDA Symp , with instantaneous spectral and bispectral features from point-process models for heartbeat dynamics. The choice of a point process framework as a physiologically inspired probabilistic method lies in its uniqueness in parameterizing heartbeat  www.nature.com/scientificreports www.nature.com/scientificreports/ dynamics continuously with any time resolution, without the need of an interpolation method. The quadratic autoregressive nonlinearities allows for the estimation of instantaneous bispectral measures and their decomposed high frequency component, which refers to the superimposed phasic activity. Time-varying information of such instantaneous phasic dynamics was quantified through central tendency (median) and AUC quantifiers. Previous studies accounted for first-and second-order moments exclusively to characterize instantaneous ANS activity 7,15,20,61-64 , therefore neglecting effective information from the superimposed phasic dynamics.
We first evaluated the proposed methodology during a paradigmatic ANS elicitation through prolonged cold-pressor test (CPT) in twenty six healthy volunteers. The significant decrease in μ RR following the CPT onset is consistent with previous studies (see 29 and references therein), and RMSSD and pNN50 also showed a significant decrease in specific time windows. Note that significant p-values were obtained from standard EDA measurements (F 4 ) at all time intervals from the CPT onset, and a delayed response in EDA Symp of about 30s from the CPT onset was also observed. We demonstrated that the sympatho-vagal balance can be better characterized by combining spectral and bispectral features derived from HRV with measures computed from EDA spectral analysis. The discriminative power of these measures is proven by the significant p-values associated with all time windows following the CPT onset, whereas the traditional LF/HF ratio was not able to discern CPT dynamics with respect to a previous resting state (Table 2 and Fig. 5). www.nature.com/scientificreports www.nature.com/scientificreports/     www.nature.com/scientificreports www.nature.com/scientificreports/ Considering the 90s time windows from the CPT onset, while the HH did not show significant changes, the HH ph (AUC) was associated with a significant p-value. This result confirms the main hypothesis of this study for which relevant information on the ANS activity can be retrieved from the superimposed phasic behavior of HRV time-varying bispectral measures. Furthermore, the proposed new indices characterizing the cardiovascular control through EDA and HRV seem to provide a more effective indicator of the sympathovagal balance then traditional indices from HRV series only 7,15,29 .
The statistical comparison between the two valence levels showed no significant differences in the proposed feature sets. However, such features allowed to discern between pleasant and unpleasant affective elicitation through an automatic decision support system. Indeed, features associated with a non-significant p-value may provide meaningful information in a machine learning framework. By applying the SVM-RFE algorithm, the discriminant power of the high oscillations of nonlinear bispectral estimates are confirmed with respect to measures obtained from first-order moments only. A balanced accuracy of 73.08% is achieved for a valence recognition through a proper selection of both spectral and higher-order spectral features and their integration with power spectral estimates of EDA and related phasic activations ( Table 4). As a comparison analysis, an HRV feature set including indices defined in the time and frequency domains (F 1 ) and related bispectral estimates (F 2 ) was fed to the classifier and an accuracy of 63.43% was obtained (Table 5), with 64% and 62.96% as positive predictive value and negative predictive value, respectively. On the other hand, a 68.48% accuracy with positive and negative predictive values of 65.52% and 71.43%, respectively, was associated with the use of a standard EDA feature set F 4 (see Table 6). By comparing results in Tables 1 and 4, we may conclude that the proposed methodology could be a viable approach to integrate HRV and EDA dynamics in case of thermal and emotional elicitations tasks.
The different number of males and females in the dataset is a limitation of this study. Moreover, although the VAS phsychological test scores have proven the strong impact of stimuli in all subjects, there is no gold-standard procedure to evoke a sympathetic elicitation both at a cardiac and at a skin sympathetic outflow level. Furthermore, the emotional responses may not be consistent across the individuals, this increasing the inter-subject variability. However, it is important to remark that we did not aim to develop the most performant emotion recognition system, but rather we aimed to devise a computational method that effectively combines HRV and EDA information to asses autonomic changes during CPT and emotional processing. We performed the cold pressor test to induce changes in the autonomic control of cardiovascular dynamics different than an emotional elicitation task. Indeed, we could have performed other autonomic maneuvers including active and passive postural changes, lower body negative pressure, valsalva maneuver, handgrip test, and others. Moreover, the administration of the cold pressor test may have caused not negligible pain that might have affected the ANS estimates derived from HRV and EDA for further analyses. We are aware that specific lab conditions are required to study complex physiological phenomena related to emotional processing, and a direct translation of the proposed methodology to emotion recognition in ecological conditions may be challenging. Another limitation of this study is that in spite of the fact that the bispectral measure HH can be considered as a reliable index of cardiovascular vagal activity, the integration of LH and LL in the defined frequency bands may contain oscillations originating in the low frequency band which may not be linked to a a sympathetic activity exclusively.
Due to the asymmetric nature of the EDA signals, asymmetrically shaped Daubechies (dbN) wavelets were widely used to analyze the raw waveforms 65,66 . In our study, db5 (chosen among db1 to db10) was the most appropriate choice of the mother wavelet for a better decomposition of tonic and phasic activities. However, a procedure to select a specific mother wavelet should be carried out in the future. Moreover, although previous studies found a correlation between EDAsymp and sympathetic activity 7,15 , sympathetic sudomotor nerve activity may be different from the cardiac sympathetic dynamics. Further studies are therefore needed to quantitatively investigate the functional relationship between different sympathetic markers from physiological series.
Despite the mentioned limitations, this study opens new horizons for the study of autonomic control mechanisms and the functional assessment of emotional responses of human. Future work will be directed to the evaluation of the proposed methodology on data gathered during orthostatisc stressors such as head-up tilt, handgrip, and lower body negative pressure, as well as from different patient population with cardiovascular or mental disorders, or autonomic dysfunctions.

Data availability
The informed consent forms signed by the subjects prevent data from being publicly available. Data may be requested via email by researchers, upon reasonable request and verification of all ethical aspects, at: gaetano. valenza@unipi.it.