Comparative evaluation of methods for the detection of electrodermal responses to multilevel intensity thermal noxious stimuli

In this study, three different methods to identify event-related electrodermal responses (ER.EDRs) were applied, comparing noxious heat stimulation at four different levels. Although the level of heat leading to subjective pain is individual, temperatures above 40 to 45 °C are normally experienced as painful due to the stimulation of heat-sensitive nociceptors. Skin conductance (SC) data from 96 healthy participants aged between 18 and 65 years old were analyzed by means of continuous decomposition analysis (CDA), discrete decomposition analysis (DDA), and trough-to-peak analysis (TTP). Several EDR features estimated from these methods were statistically compared. Within the analysis window following each single stimulation, CDA and DDA methods identified more electrodermal responses between 1 and 9 s after noxious stimulation than TTP for all stimulation intensities. However, the occurrence rates of at least one ER.EDR after noxious stimulation were similar among the three methods and tended to increase with increasing intensities. Among the common features, amplitude sum had better discriminative power for differentiating noxious stimulation intensities regardless of the method. The results suggest that all investigated methods performed similarly in identifying electrodermal changes in response to high-intensity thermal noxious stimuli. In the experimental conditions of cutaneous heat stimulation of nociceptors, good discrimination among stimulation intensities was found using the amplitude sum feature. In conclusion, although CDA and DDA were more sensitive than TTP for identifying ER.EDRs, neither CDA nor DDA brought further discriminative power to the study of noxious stimulation intensities when compared with the traditional TTP method.


Introduction
Recordings of electrodermal response (EDR) are mainly based on the measurement of skin conductance. Skin conductance changes with the activity of the sweat glands. This response of the skin is referred to as electrodermal activity (EDA), and it is considered to reflect the activity of the sympathetic branch of the autonomic nervous system. EDR is a way to measure sympathetic responses to mental or physical stimulation, particularly in response to noxious stimulation (Bach 2014;Leknes et al. 2008). EDR can be described as a slowly varying monophasic signal that may change to a fast increase in its amplitude with slow recovering tails. These abrupt deflections can be related to either internal arousal processes or external stimuli (event-related electrodermal responses (ER.EDRs)). Due to the slow recovery process following an EDR, studies using varying interstimulus intervals (ISIs) showed that responses following stimulation with ISIs as low as 2 s cannot be properly handled by the traditional trough-to-peak analysis (TTP) because such stimulation elicits responses that are superimposed on one another and appears as a single EDR in the raw signal (Benedek and Kaernbach 2010a, b). To overcome this problem, decomposition methods by means of deconvolution can be applied to sudomotor nerve activity (SNA) data (Bach et al. 2010;Bach and Friston 2013;Benedek and Kaernbach 2010b, a).
SNA and EDR may serve as important sources of information, especially in scenarios involving non-communicative patients in which pain is often inadequately identified and treated (Lichtner et al. 2014;Lucas et al. 2016;Walther-Larsen et al. 2017). Experimental noxious stimulation of healthy participants is used to better understand the nociceptive mechanisms involved in pain perception. This stimulation can be done by means of different agents, including intense pressure, subcutaneous injection of chemical substances, electrical stimulation, and extreme temperatures in both cold and hot domains. One of the most used modalities is the heat stimulation, in which temperatures above 40 to 45°C are expected to activate thermal nociceptors and lead to the perception of pain depending on the frequency of action potentials, temporal summation, and central modulation of afferent signals (Dubin and Patapoutian 2010).
Although some studies have mentioned the low specificity and low predictive values of EDR analysis in assessing pain phenomena or the response to noxious stimulation (Ahmed et al. 2015;Ledowski et al. 2009;Strehle and Gray 2013;Valkenburg et al. 2012), many studies have analyzed EDRs by themselves or in combination with other physiological responses to monitor the pain sensation and try to provide more objective measurements (Ahmed et al. 2015;Günther et al. 2016;de Jesus et al. 2015;Kächele et al., 2015a, b;Loggia et al. 2011;Sabourdin et al. 2013;Storm 2013;Strehle and Gray 2013;Treister et al. 2012;Walter et al. 2014;Werner et al. 2014).
In this research, TTP and decomposition analyses are applied to signals resulting from electrodermal activity with the aim of detecting and discriminating noxious thermal stimuli at several intensities.

Database and participants
This work was developed using previously recorded data from the Biovid Pain Database project (ethics committee approval: 196/10-UBB/bal, Ulm University, Ulm, BW, Germany) (Walter et al. 2013).
This database consists of a multimodal dataset obtained by monitoring several biopotentials during the application of noxious thermal stimulation on the forearm skin at four different intensities. For this study, only electrodermal signals were investigated. Data from 96 healthy participants equally distributed by gender and aged between 18 to 65 years old (40.9 ± 14.9 years, mean ± standard deviation) were used.

Noxious stimulation
Thermal stimuli were applied to the backside of the right forearm while the participant sat in a chair with both arms resting on a table. The experiment was performed using a thermode (PAHTWAY, Medoc, Israel). For each participant, the subjective heat pain threshold (P1) and tolerance limit (P4) stimulus intensities were obtained by continuously increasing the thermode's surface temperature. P1 is defined as the lowest stimulus intensity able to cause pain, while P4 is defined as the maximum temperature the participant can stand. A maximum temperature of 50.5°C was not exceeded in order not to harm skin integrity. Thereafter, two equally distributed intermediate temperatures (P2 and P3) were calculated in order to obtain four distinct stimulation intensities. Finally, each temperature (P1, P2, P3, and P4) was randomly applied 20 times and sustained for 4 s, with randomized pauses of 8 to 12 s between consecutive stimuli (baseline temperature of 32°C) (Walter et al. 2013).

Data acquisition
EDR was recorded using a Nexus-32 amplifier (Humakarigar Pvt. Ltd., India), with a sampling frequency of 512 Hz, and event-related data (i.e., thermal stimulus) were simultaneously recorded using Biotrace+ software (Mind Media, Netherlands).

Data analysis
ER.EDRs were analyzed in MATLAB R2015a (MathWorks Inc., USA) using the open-source software Ledalab v3.4.8 (available at http://www.ledalab.de), with 8-s acquisition window (between 1 and 9 s; see Eq. 3). Data were filtered with 1 st order Butterworth low-pass filter with a 5 Hz cutoff frequency, downsampled to 64 Hz, and analyzed using the three methods described below. A block diagram depicting the data analysis sequence is provided in Fig. 1.

Trough-to-peak analysis
TTP analysis is the traditional and simplest technique for identifying EDRs. It consists of running peak detection algorithms over the EDA signal to detect positive deflections on it. Thresholds of 0.03 μS or 0.05 μS are normally used (Boucsein 2012;Boucsein et al. 2012) to detect significant EDRs. Commonly estimated features for the characterization of EDRs include amplitude, rise time, half recovery time, and latency Lukander 2016).

Discrete decomposition analysis
Decomposition techniques developed by Benedek and Kaernbach (2010b) are based on the works of Alexander et al. (2005) and Lim et al. (1997), in which the electrodermal signal is composed of a slowly varying tonic activity and a phasic component, which is expressed as the convolution of the underlying SNA and an impulse-response function (IRF). The main advantage of discrete decomposition analysis (DDA) is the possibility to detect deviations in the EDR shape with regard to a fixed response function (Benedek and Kaernbach 2010b).
Therefore, the authors assume that the electrodermal signal consists of a non-negative driver characterized by compact non-overlapping positive deflections over a zero-valued baseline that represents the SNA, an IRF, while the remainder that accounts for variation in EDR shapes is supposed to be related to the pore opening processes of the sweat glands (Benedek and Kaernbach 2010b). The IRF-b(t)-consists of a biexponential function with two time constants, τ 1 and τ 2 , which account for a steep onset and a slow recovery: A full description of the method can be found elsewhere (Benedek and Kaernbach 2010b).

Continuous decomposition analysis
Continuous decomposition analysis (CDA) was developed by Benedeck and Kaernbach and is fully described in their work (Benedek and Kaernbach 2010a). The EDA signal is estimated as follows: in which SCL (skin conductance level) is the tonic component of the EDA signal, and the IRF is given in (1).

Description of the set of features
The following steps summarize how features are calculated: & Set the beginning of the analysis window to 1 s after the beginning of stimulation. & Set the end of the analysis window to 9 s after the beginning of stimulation. & Find all local minima followed by local maxima within the analysis window, if any. & Exclude pairs of local minima and local maxima that have differences less than 0.01 μS.
& Proceed with analysis for any pair of local minima and local maxima with differences greater than 0.01 μS.
Let S be the windowed signal as a function of the finite time set T, which are both ordered sets, described as follows: Given π and υ, V and P are ordered sets containing valleys and peaks of S, and TV and TP are ordered sets containing the timestamps in which valleys and peaks occur, so that: Whenever V and P are non-empty sets, significant responses can be retrieved by computing the relative amplitudes between each identified valley and the following peak whose differences are larger than a threshold, θ: given: After detecting significant peaks in the windowed signal, a number of features described in Table 1 are computed. Only EDRs with latencies between 1 and 9 s and a relative amplitude larger than 0.01 μS were considered ER.EDRs. The occurrence rate (OCR) refers to the percentage of cases in which at least one ER.EDR is detected in a specific condition, so that if there were 20 stimuli at P1 level and significant ER.EDRs were found in only five of them, then OCR = 25%. This feature indicates the ability of the analysis to detect ER.EDR due to noxious stimulation.
The following features are automatically calculated by the Ledalab software: the number of ER.EDRs or number of skin conductance responses (NSCR), Latency, and amplitude sum (AmpSum). These features were computed for all three methods (CDA, DDA, and TTP). The sum of identified ER.EDRs (AreaSum) was only computed for DDA, and average phasic driver (SCR), area of phasic driver (ISCR), and maximum value of detected ER.EDRs (PhasicMax) were only computed for CDA.

Statistical analysis
Statistical analyses were performed using R Software version 3.5.2 (R Core Team 2017). Since Lilliefors normality test rejected the hypothesis of normal distribution of all variables, data from NSCR, Latency, Ampsum, and OCR features were transformed using the aligned rank transformation in order to enable the study of interactions in 4 (intensity: P1 vs. P2 vs. P3 vs. P4) by 3 (method: CDA vs. DDA vs. TTP) betweensubject design. One-way and two-way analyses of variance were used on the aligned ranks to assess the main effects in which R is the set of the relative amplitudes of the responses within the analysis window using either TTP, DDA, or CDA.
AreaSum dda (μS*s) Sum of the identified ER.EDR areas using DDA in which S is the convolved signal of the phasic driver and the IRF, and t 1 and t n are the bounds of the analysis window. SCR cda (μS) Average phasic driver estimated using CDA in which s i is the ith sample of the phasic driver signal. ISCR cda (μS*s) Estimated area of the phasic driver using CDA in which s x are elements of the phasic driver signal (S), and t 1 and t n are the bounds of the analysis window.
CDA, continuous decomposition analysis; DDA, discrete decomposition analysis; ER.EDRs, event-related electrodermal responses; TTP, trough-to-peak analysis Fig. 1 Block diagram depicting the sequence of steps employed for data analysis: CDA, continuous decomposition analysis; DDA, discrete decomposition analysis; EDA, electrodermal activity; TTP, trough-to-peak analysis and interactions, respectively, as suggested in the literature (Wobbrock et al. 2011). The other features were not transformed, and the Friedman test was used to assess the main effects of stimulus intensity. Bonferroni correction for posthoc test was used to discriminate significant differences between methods, intensities, and their interactions. The significance level was set at 0.05 for all comparisons.

Results
Data from one participant were excluded due to a corrupted EDA signal. The final dataset consisted of 95 complete EDA signals.
The NSCR feature presented significant main effects of both intensity, F(3.23) = 314.94, p < 0.0001, and method, F(2.23) = 276.12, p < 0.0001. However, these main effects were qualified by a significant interaction between the two factors, F(6.23) = 4.26, p = 0.0003. Bonferroni-adjusted comparisons indicated that for the same intensity, NSCR values obtained with CDA and DDA were larger than those using TTP for all intensities at a significance level lower than 0.01%. NSCR values obtained with CDA and DDA only differed for intensities P1 (p = 0.0051) and P2 (p = 0.0067). For the same method, NSCR values obtained with DDA at P3 were larger than those at P1 (p = 0.0261) and at P2 (p = 0.0446). However, no significant differences for NSCR values obtained with CDA or TTP methods were observed. The probability density estimates of NSCR for each method of analysis and stimulus intensity are shown in Fig. 2.
The OCR feature of significant responses presented significant main effects of the intensity factor, F(3.11) = 41.31, p < 0.0001, in which P3 and P4 values were larger than P1 (p < 0.0001 and p < 0.0001, respectively) and P2 (p = 0.0110 and p < 0.0001, respectively) values, but demonstrated no main effects of the method, F(2.11) = 0.89, p = 0.4101. Furthermore, no significant interactions were found between the two factors, F(6.11) = 1.61, p = 0.1414. The probability density estimates of OCR for each method of analysis and stimulus intensity are shown in Fig. 3.
The AmpSum feature presented significant main effects of both intensity, F(3.23) = 516.52, p < 0.0001, and method, F(2.23) = 6.28, p = 0.0019. However, these main effects were qualified by a significant interaction between the two factors, F(6.23) = 8.07, p < 0.0001. Bonferroni-adjusted comparisons indicated that for the same method, all intensities were significantly different for all three methods, with larger values for the higher intensities at a significance level equal to or less than 0.01%. For the same intensity, AmpSum values obtained with CDA and DDA were larger than those using TTP for Fig. 2 Probability density estimates of the number of skin conductance responses (NSCR) and their medians (vertical lines) after applying 80 random stimuli with four increasing intensities (P1, P2, P3, and P4, 20 stimuli each) detected by using continuous decomposition analysis (CDA), discrete decomposition analysis (DDA), and trough-topeak analysis (TTP) intensities P1, P2, and P3 at a significance level equal to or lower than 0.01%, but their values were lower than those using TTP for P4 intensity (p < 0.0001). Furthermore, no significant differences between AmpSum values obtained with CDA and DDA for any of the four intensities were observed; however, in P1 (p = 0.0216), the values obtained with CDA were lower than those using DDA. The probability density estimates of AmpSum for each method of analysis and stimulus intensity are shown in Fig. 4.
The Latency feature presented significant main effects of both intensity, F(3.13) = 69.21, p < 0.0001, and method, F(2.13) = 90.91, p < 0.0001. However, these main effects were qualified by a significant interaction between the two factors, F(6.13) = 5.31, p < 0.0001. Bonferroni-adjusted comparisons indicated that at P4 intensity, Latency values obtained with CDA and DDA were lower than those using TTP (p = 0.0409 and p = 0.0068). For the same method, differences between intensities were only significant for Latency values obtained with TTP, in which P1 and P2 values were lower than P4 values (p = 0.0078 and p = 0.0074). The probability density estimates of Latency for each method of analysis and stimulus intensity are shown in Fig. 5.

Discussion
In this research, TTP and decomposition analyses were applied to signals resulting from electrodermal activity. The aim was to detect and discriminate noxious thermal stimuli at several intensities.
Although CDA and DDA presented significant improvements in the count of event-related EDRs when compared with TTP for all stimulus intensities, this fact did not present significant improvements in the OCR of ER.EDRs regardless of the intensity. This finding means that more EDRs were identified, but solely where at least one single response had already been identified using TTP analysis. Fig. 3 Probability density estimates of occurrence rates (OCR) of ER.EDRs and their medians (vertical lines) after applying 80 random stimuli with four increasing intensities (P1, P2, P3, and P4, 20 stimuli each) detected by using continuous decomposition analysis (CDA), discrete decomposition analysis (DDA), and trough-to-peak analysis (TTP) The second point was to assess whether CDA and DDA could provide better discrimination among noxious stimulation intensities than TTP. Our results showed an improvement in the number of identified responses, as well as in the amplitude features. However, these results did not support improvements in discrimination power among the four intensities when compared with TTP.
Considering the NSCR feature, which the literature usually considers a feature relevant to pain evaluation, it tends to increase with increased pain or increased stimulation intensity. Clinical studies have found from a moderate correlation with numeric rating scales in postoperative inpatients (Rago et al. 2015) to no correlation with behavioral scales in neonates during the heel prick procedure (de Jesus et al. 2015), while an experimental study using thermal stimulation found moderate discrimination power among different stimulation intensities (Treister et al. 2012) in healthy participants. Similarly, a recent study showed that an increased NSCR is more closely related to the electric noxious stimulation than to positive or negative emotions elicited by pictures (Günther et al. 2016). Our findings corroborate these results, and, although values obtained with CDA and DDA methods were significantly larger than those obtained with TTP method for all stimulus intensities, if one intends to discriminate intensities of noxious stimuli by means of NSCR values, DDA method performed slightly better than CDA and TTP methods.
Most previous studies corroborate the idea that amplitude values also tend to increase with increased thermal stimulation intensities but without good decimation power related to stimulus intensity or severity of pain (Loggia et al. 2011;Treister et al. 2012;de Jesus et al. 2015). In contrast, the study of Walter et al. (2014) found a set of new EDR features related to amplitude that could provide discriminative power among four heat noxious stimulation intensities. In the present work, only the AmpSum feature could accurately discriminate the four noxious stimulation intensities, regardless of the method used to compute it.
Other measures that are somehow related to the signal's amplitude (PhasicMax c d a , SCR c d a , ISCR c d a , and AreaSum dda ) and were calculated only for the decomposition methods also increased with increased noxious stimulation intensities. Although they performed better than other common features evaluated in this paper, they were only able to discriminate the higher stimulation intensities from the lower ones, but none of them could discriminate all four intensities of noxious stimulation.
Finally, CDA and DDA presented smaller Latency values than TTP, which may account for the identification and Fig. 4 Probability density estimates of amplitude sum (AmpSum) of the ER.EDRs and their medians (vertical lines) after applying 80 random stimuli with four increasing intensities (P1, P2, P3, and P4, 20 stimuli each) detected by using continuous decomposition analysis (CDA), discrete decomposition analysis (DDA), and trough-to-peak analysis (TTP) discrimination of hidden or superposed responses at the beginning of the stimulation phase. Although a subtle difference has been observed among the methods, only TTP could discriminate the stimulus intensities. In addition, no recent studies were found in which a relationship between latency and stimulus intensity has been investigated.
Considering these recent papers, it may be noted that amplitude values and NSCR can provide some reliable information on nociceptive perception. Using sustained thermal stimulation for a few seconds (approximately 6 s) and a comprehensive response window, our results partially support this statement. However, one drawback of using ER.EDRs to study how the human body reacts to the noxious heat stimulation is that it still relies on the lack of specificity and on the low predictive power for a clinical application to quantify how strongly patients respond to different modalities of noxious stimuli. Even when using new methods to estimate the SNA instead of using EDR data itself, no advances were achieved in terms of discrimination of noxious stimulation intensity. All three methods presented similar discrimination power for the common features among them. Some authors (Loggia et al. 2011;Walter et al. 2014) suggest that searching for new features may provide better insight into the study of sympathetic responses to noxious stimulation, responses which may be related to the subjective experience of pain. In our study, the four features that were only computed for the CDA or DDA methods presented some tendency but poor discrimination power among the four stimulus intensities. In contrast, only the common feature AmpSum could discriminate all four intensities of noxious stimulation, regardless of the method. Therefore, future works should keep looking for new significant features while using more robust signal processing methods, such as CDA and DDA, but should not be restricted to them.
Regarding decomposition methods, to the best of our knowledge, at least two ready-to-use tools could be further analyzed to check our results on decomposition methods. Those include the PsPM, which uses a canonical skin conductance response function and a general linear convolution model (Bach et al. 2010), and the CvxEDA, which uses Bayesian statistics, a convex optimization approach, and an IRF based on an infinite impulse response (IIR) (Greco et al. 2016). Although the developers of these other methods have reported significant improvements in EDR analysis when compared with Ledalab software, we opted to use the Ledalab software, which combines the three different analytic methods, for our first trial. Therefore, it can be emphasized that different methodologies should be further evaluated, and generalizations to other decomposition methods should be carefully done or altogether avoided. Furthermore, with regard to the differences among the three methods, using the specified parameters, the CDA and TTP methods were faster than DDA.  5 Probability density estimates of the latency of the ER.EDRs and their medians (vertical lines) after applying 80 random stimuli with four increasing intensities (P1, P2, P3, and P4, 20 stimuli each) detected by using continuous decomposition analysis (CDA), discrete decomposition analysis (DDA), and trough-to-peak analysis (TTP)

Conclusion
In summary, it can be concluded that all the investigated methods can be equally suitable to identify electrodermal skin responses to high-intensity noxious stimuli. CDA and DDA were more sensitive at identifying EDRs compared with the TTP method, indicating that they more efficiently identify slight activations of the sudomotor nerves that are not obvious or easily identified in the raw EDR signal. Although EDR occurrence rates tended to increase with increased stimulus intensity regardless of the method of analysis, ranging from 40 to 78% of the time, no differences among them for the same intensities were observed. The amplitude features of ER.EDRs using any of the evaluated methods can provide better discriminative power than comparing the number of evoked responses among the four stimulus intensities proposed here.