Elementary integrate-and-fire process underlies pulse amplitudes in Electrodermal activity

Abstract

Electrodermal activity (EDA) is a direct read-out of sweat-induced changes in the skin’s electrical conductance. Sympathetically-mediated pulsatile changes in skin sweat measured as EDA resemble an integrate-and-fire process, which yields an inverse Gaussian model as the inter-pulse interval distribution. We have previously showed that the inter-pulse intervals in EDA follow an inverse Gaussian distribution. However, the statistical structure of EDA pulse amplitudes has not yet been characterized based on the physiology. Expanding upon the integrate-and-fire nature of sweat glands, we hypothesized that the amplitude of an EDA pulse is proportional to the excess volume of sweat produced compared to what is required to just reach the surface of the skin. We modeled this as the difference of two inverse Gaussian models for each pulse, one which represents the time required to produce just enough sweat to rise to the surface of the skin and one which represents the time requires to produce the actual volume of sweat. We proposed and tested a series of four simplifications of our hypothesis, ranging from a single difference of inverse Gaussians to a single simple inverse Gaussian. We also tested four additional models for comparison, including the lognormal and gamma distributions. All models were tested on EDA data from two subject cohorts, 11 healthy volunteers during 1 hour of quiet wakefulness and a different set of 11 healthy volunteers during approximately 3 hours of controlled propofol sedation. All four models which represent simplifications of our hypothesis outperformed other models across all 22 subjects, as measured by Akaike’s Information Criterion (AIC), as well as mean and maximum distance from the diagonal on a quantile-quantile plot. Our broader model set of four simplifications offered a useful framework to enhance further statistical descriptions of EDA pulse amplitudes. Some of the simplifications prioritize fit near the mode of the distribution, while others prioritize fit near the tail. With this new insight, we can summarize the physiologically-relevant amplitude information in EDA with at most four parameters. Our findings establish that physiologically based probability models provide parsimonious and accurate description of temporal and amplitude characteristics in EDA.

Introduction

Sweat gland activity is used to assess sympathetic nervous system activity in applications such as lie detector tests and neuromarketing [1]. Sympathetic activation is also known as the “fight or flight response”, which is induced by states such as stress, anxiety, and pain [1]. Electrodermal activity (EDA) measures the second-to-second electrical conductance of the skin to capture sweat gland activity. As stimulation of sweat glands increases due to stress or pain for example, more sweat is produced, which increases the electrical conductance of the skin. EDA is typically divided into two components [1]. The first is a baseline or tonic component which drifts gradually over minutes and is thought to represent ambient conditions which contribute to baseline level of filling of the glands. The second is the phasic component, which rides on top of the tonic and consists of pulsatile sweat release events. These pulsatile sweat release events have a timescale of a few seconds and are thought to correspond more closely to sympathetic nervous system activity [1]. There is growing interest in the development of algorithms to accurately characterize changes in emotional and physiologic states from EDA.

Recent EDA modeling efforts fall into several categories, as outlined in [2]. Increasingly advanced and diverse tools from signal processing domains are being employed to design new decomposition methods to separate EDA into tonic and phasic components [3–11]. However, each of these methods yields different results on the same datasets [2] and none have physiological validation mechanisms. Some approaches involve designing frequency domain measures to analyze EDA, based on analogous methods for frequency domain analysis of heart rate variability [12–15]. In the subset of approaches in which pulse amplitude is examined specifically, the pulse amplitude is assumed to correlate with stimulus intensity in controlled experimental settings [3–11,16–20]. Where state space approaches are used to model both pulse occurrence and amplitude, the amplitudes are assumed to follow a Gaussian distribution [21,22].

Our previous analyses have showed that the inter-pulse interval distribution in EDA data follows an inverse Gaussian distribution, which agrees with a model of the rise of sweat through the gland to the skin surface as an integrate-and-fire process, specifically a Gaussian random walk with drift diffusion [23,24]. We showed that deviations from the inverse Gaussian due to recording across many sweat glands tend toward right-skewed heavier tailed distributions, such as the lognormal [23]. Using these insights, we further defined a low-order paradigm for verifying the physiologic structure in EDA that includes a framework for goodness-of-fit analysis [25,26].

However, temporal information is not the only information in phasic EDA. Each pulse occurs not only at a discrete point in time, but also with a specific amplitude [1]. Existing algorithms for phasic EDA analysis assume that the amplitude of each pulse has a one-to-one association with the intensity of the stimulus driving it [3–11,16–20]; however, previous physiologic experiments have shown that the background level of nerve activity and baseline level of filling of the sweat glands can alter pulse amplitude even in the face of an unchanging stimulus intensity [27–29]. In this work, we propose a model for pulse amplitudes that uses the same insight about integrate-and-fire physiology of sweat glands as for temporal information. We hypothesize that the amount of sweat produced in a pulse relates directly to stimulus amplitude. However, we postulate that the observed amplitude of the pulse relates to how much more sweat is produced than what is required to reach the surface of the skin, which also accounts for the role played by the background filling level of the sweat glands. Therefore, we model the amplitude of a pulse as the difference of actual amount of sweat produced and the baseline filling level.

We implement four different simplifications of our model in two different subject cohorts. The four simplifications tested were a simple inverse Gaussian, a three-parameter inverse Gaussian with the third parameter serving as a location shift, a single difference of an inverse Gaussian and a Gaussian, and a single difference of two inverse Gaussians. We show, using a goodness-of-fit analysis, that each simplification balances the important characteristics of the model differently. The simple inverse Gaussian model fits the mode of the pulse amplitude distribution well, while the difference models better capture the tail. The three-parameter model seems to balance both.

Important advances we report are a set of low-order physiology-based point process models for pulse amplitudes in phasic EDA that work synergistically with our existing models for temporal information. Using both together, we can extract all relevant information from phasic EDA in statistically rigorous way. The balance of this paper is organized as follows. In Materials & Methods, we derive our hypothesis about pulse amplitudes from the integrate-and-fire physiology, outline four statistical models to capture it that involve the inverse Gaussian as well as four alternatives for comparison. In Results, we use these models in the analysis of EDA pulse amplitudes recorded from 22 subjects across two different subject cohorts, one while awake and at rest and the other under controlled propofol sedation. The Discussion describes the implications of our findings for future basic science and translational studies.

Materials & methods

Anatomy and physiology

We review the anatomy and physiology of sweat production in the skin in detail in the S2 Appendix. The pulsatile changes in conductance measured in the skin are referred to as ‘pulses’ in this paper. Existing algorithms for EDA analysis typically assume that pulse amplitude can be explained solely by the intensity of the stimulus, and therefore they rely on the pulse amplitude to directly infer stimulus amplitude [3–11,16–20]. However, physiologic experiments done by Wallin et al. in the 1990s demonstrated that modulating the background alone could result in pulses of varying amplitudes, even if stimulus intensity was held constant [27–29]. Therefore, interpreting pulse amplitude in light of stimulus intensity alone, especially in a context in which background cannot be held constant, can be misleading. Any physiologically viable model for pulse amplitude must also account for the background.

Building upon the integrate-and-fire model we postulated for temporal information in sweat gland bursts, we hypothesized that the amplitude of a pulse is proportional to the excess volume of sweat produced compared to what would be required to reach the surface of the skin, where the total volume of sweat produced relates to the intensity of the stimulus, and the amount of sweat required to reach the surface of the skin is determined by the background filling level. The background filling level is affected by tonic EDA, spontaneous activity, reabsorption rate in the duct of the sweat gland, each individual’s autonomic ‘excitability’, and the conductance properties of their skin [27–32]. It also varies from sweat gland to sweat gland.

Statistical model

By assuming a relatively linear relationship between the volume of sweat in the sweat glands and measured electrical conductance across the skin, and a relatively constant rate of sweat production once stimulated, we can relate measured electrical conductance across the skin to the times taken to secrete the required volume of sweat. Both assumptions are simplifications to the true microfluidic properties, made across the aggregate of hundreds of sweat glands [30]. With the help of these assumptions, we hypothesize that the amplitude of each pulse can be modeled as the difference of two processes: one integrate-and-fire process to reach the surface of the skin (the minimum amount of sweat production required), and a second integrate-and-fire process to reach up to a portion of the maximal capacity of the gland (the actual sweat production resulting from the intensity of the stimulus). However, each pulse may have a different stimulus intensity and background filling level, so the two processes are not identical across pulses. Since they are both integrate-and-fire processes, we hypothesize that each pulse amplitude can be modeled as the difference of two inverse Gaussian distributions as shown in Fig 1A [33–35].

[Figure omitted. See PDF.]

Fig 1. Schematic of our physiological hypothesis and statistical model simplifications.

(A) We hypothesize that the amplitude of a pulse is related to the excess volume of sweat produced relative to what is required to reach the surface of the skin, which translates to the difference of first passage times between two integrate-and-fire processes, yielding the difference of two inverse Gaussian distributions. (B) We arrived at a series of simplifications of our model based on the statistical properties observed. IG = inverse Gaussian.

https://doi.org/10.1371/journal.pcbi.1009099.g001

Since fitting this model for each pulse individually is an under-constrained problem, we proposed four viable simplifications across a single subject’s dataset (Fig 1B):

1. a single difference of inverse Gaussians (IG-IG),

2. a single difference between an inverse Gaussian and a Gaussian (IG-G),

3. a single 3-parameter inverse Gaussian with the third parameter being a location shift (3IG), and

4. a single simple inverse Gaussian model (SIG).

The first simplification is the most obvious place to start, but preliminary results indicated that the second inverse Gaussian actually approached a very narrow Gaussian distribution [35], leading to the second and third simplifications. The fourth simplification is the simplest of all four. Our previous work indicated that the inverse Gaussian was the best integrate-and-fire model for EDA data across subjects [23], and therefore we only included the inverse Gaussian rather than the larger family of all integrate-and-fire models. We also compared other families of right-skewed models, such as the lognormal, exponential, and gamma. We performed a goodness-of-fit analysis for all models with quantile-quantile (QQ) plots and rescaled QQ plots. QQ plots show the quantiles of one distribution against another, in this case comparing the theoretical and empirical distributions [36,37]. We also calculated Akaike’s Information Criterion (AIC) and the mean and maximum distance from the 45-degree line on the QQ plots [36,37]. Models 1 and 2 implicitly assumed independence of the two distributions involved; the validity of this assumption was left to verification by the goodness-of-fit analysis.

Experimental data

We tested all of the models on two subject cohorts, collected at different times using different equipment and under different conditions. The first cohort of data is EDA data we previously collected from 12 healthy volunteers between the ages of 22 and 34 (6 males) while awake and at rest [23]. The study was approved by the Massachusetts Institute of Technology (MIT) Committee on the Use of Humans as Experimental Subjects. All subjects provided written informed consent. Approximately one hour of EDA data was collected at 256 Hz from electrodes connected to the second and fourth digits of each subject’s non-dominant hand using the Thought Technology Neurofeedback System [38]. Data from the electrodes were fed into an encoder connected to a laptop in the neighboring room via fiber optic cable. We monitored the data as it was collected for the entire hour. Subjects were seated upright and instructed to remain awake. They were allowed to read, meditate, or watch something on a laptop or tablet, but not to write with the instrumented hand. One subject’s data were not included in the analysis because we learned after completing the experiment that the subject occasionally experienced a Raynaud’s type phenomenon, which would affect the quality of the EDA data. Data from the remaining 11 subjects were analyzed.

The second data cohort consists of EDA recorded from eleven healthy volunteers during a study of propofol-induced unconsciousness [39]. The protocol was approved by the Massachusetts General Hospital (MGH) Human Research Committee. All subjects provided written informed consent. For all subjects, approximately 3 hours of data were recorded while using target-controlled infusion protocol. The data collection, including experimental setup, is described in detail in [39]. EDA data were recorded using the BedMaster system by placing 2–3 electrodes on left hand of each subject [40]. The infusion rate was increased and then decreased in a total of ten stages of roughly equal lengths to achieve target propofol concentrations of: 0 mg/kg/hr, 1, 2, 3, 4, 5, 3.75, 2.5, 1.25, 0. All data were analyzed using Matlab R2019a. [41].

Data preprocessing and pulse selection

Preprocessing consisted of two major steps, 1) detecting and removing artifacts and 2) isolating the phasic component. Both have been described in previously in [26]. Because of the level of high frequency noise seen in the recording equipment used for the propofol data, those data were additionally low-pass filtered with a cutoff of 3 Hz after artifact removal.

Pulse selection was done using the methodology described in [25] in which the fits of four right-skewed models were used to select the best prominence threshold at which to extract pulses. Prominence is a locally adjusted amplitude measure computed using the findpeaks algorithm in Matlab. This algorithm adjusts the amplitude of each peak in the signal as the height above the highest of neighboring “valleys” on either side. The valleys are chosen based on the lowest point in the signal between the peak and the next intersection with the signal of equal height on either side. Since the same pulse selection framework was followed on the same two cohorts of data, the pulses selected for each subject were also the same as in [25]. The final thresholds used for each subject and temporal properties of the pulses selected can be found in [25]. For this paper, we used the prominence of the extracted pulses as the measure of pulse amplitude.

Statistical model fitting and comparison

We fitted eight models to each subject’s dataset of extracted pulse amplitudes (prominences). The first four were the four simplifications of our hypothesis, and the other four were other models for comparison. The eight models fitted were:

1. a single difference of inverse Gaussians (IG-IG),

2. a single difference between an inverse Gaussian and a Gaussian (IG-G),

3. a single 3-parameter inverse Gaussian with the third parameter being a fitted location shift (3IG with fitted shift),

4. a single simple inverse Gaussian model (SIG),

5. a single lognormal model (L),

6. a single gamma model (G),

7. a single exponential model (E), and

8. a single 3-parameter inverse Gaussian with the location shift parameter set to the prominence threshold used to extract pulses (3IG with known shift).

The closed form densities for Models 3–8 are in Table 1. We fitted models 4–7 by maximum likelihood [42]. We fitted models 1–3 and 8 by method-of-moments [43], due to a lack of closed form solutions for maximum likelihood estimates of the parameters. The derivation of method of moments estimates for Models 1 and 2 is detailed in S3 Appendix. Models 1 and 2 do not have closed form densities, which we have noted in Table 1. Method of moments estimates for Model 3 (also used for Model 8) were given in [35]. Model 8 was included to verify that the process of pulse extraction using a prominence threshold did not skew the pulse amplitude results. We assessed goodness-of-fit by Akaike’s Information Criterion (AIC) and QQ plots. The AIC is defined aswhere is the likelihood evaluated at the maximum likelihood estimate of the parameters and p is the number of parameters. A lower AIC indicates a better fit. For the models fitted by method of moments, we estimated the log likelihood numerically. AIC prioritizes efficiency of the model.

[Figure omitted. See PDF.]

Table 1. Model descriptions for Models 1–8.

https://doi.org/10.1371/journal.pcbi.1009099.t001

We also plotted QQ plots and rescaled QQ plots. For rescaled QQ plots, both model and empirical quantile values were rescaled byRescaled QQ plots were used for more uniform visualization of the data across all quantiles. From the QQ plots (not rescaled), we calculated the mean and maximum perpendicular distances between the plotted fit and the 45-degree line, which represents a perfect fit. A lower mean distance indicates a better average fit across all quantiles, while a lower maximum distance indicates a better worst case fit (a better worst-fitting point).

Results

Extraction of pulses

In the awake and at rest cohort, the number of pulses extracted per subject across one hour ranged from 97 to 324 using prominence thresholds ranging from 0.0025 to 0.027. In the propofol sedation cohort, the number of pulses extracted per subject across 3–4 hours ranged from 383 to 1250 using prominence thresholds ranging from 0.02 to 0.055.

Findings from statistical model comparison

The detailed results for the simplification models (Models 1–4) are in Tables 2–3, and the detailed results for the other models (Models 5–8), which performed poorly overall, are Table A in S1 Appendix and Table B in S1 Appendix. Based on AIC, in the awake and at rest cohort (Tables 2 and A), SIG is the best model for 9 out of the 11 subjects (S3-S11), lognormal for one subject (S2), and the 3IG with known location shift (Model 8) for one subject (S1). The SIG was always the best of the four simplifications (Models 1–4). In the propofol sedation cohort (Tables 3 and B), SIG was the best model for 8 of 11 subjects (P1-P7, P11), 3IG with fitted shift for one (P8) and lognormal for two (P9, P10).

[Figure omitted. See PDF.]

Table 2. Model fit results for awake and at rest cohort for Models 1–4.

https://doi.org/10.1371/journal.pcbi.1009099.t002

[Figure omitted. See PDF.]

Table 3. Model fit results propofol sedation cohort for Models 1–4.

https://doi.org/10.1371/journal.pcbi.1009099.t003

Based on mean distance, in the awake and at rest cohort, SIG was the best model for 8 of the 11 subjects (S1, S3-S8, S11), lognormal for one subject (S2), and 3IG with fitted shift for two subjects (S9, S10). In the propofol sedation cohort, 3IG with fitted shift was the best model for 7 of the 11 subjects (P2-P5, P7, P8, P10) and SIG for the other 4 (P1, P6, P9, P11).

Based on max distance, in the awake and at rest cohort, one or more of the two difference models and the 3IG with fitted shift model were the best across 10 of the 11 subjects (more than one model performed equally well in most cases). The 3IG with fitted shift model was one of the best models for 9 of the 11 subjects (S1-S4, S6-S9, S11), the IG-G model for 5 of the 11 subjects (S3, S5-S7, S11), and the IG-IG model for 8 out of 11 subjects (S1, S3, S4, S6-S9, S11). The exponential was the best model for the remaining subject (S10). In the propofol sedation cohort, the 3IG with fitted shift was the best for 7 out of 11 subjects (P2, P5, P7-P11) and the IG-G for the remaining 4 (P1, P3, P4, P6).

Overall, all simplifications (Models 1–4) were reasonable models for the data and outperformed other distributions (Models 5–8). Each of the simplifications of our hypothesis (Models 1–4) prioritized different aspects of the fit (Figs 2–5 and Figs A-AN in S1 Appendix). The maximum distance from the 45-degree line seemed to occur at high quantiles, which was the right tail for most models. The fit in this region was prioritized by the difference models, Models 1 and 2 (IG-IG and IG-G). Models 1 and 2 fitted the right tail of the distribution by sacrificing some of the fit near the mode of the distribution, reflected in mean distance and AIC, especially as compared to SIG (Model 4). The 3IG with fitted shift (Model 3) seemed to balance the fit of both mode and tail of distribution reasonably. These results suggest that SIG was the best model in terms of efficiency, but if overall quality of fit was prioritized, the 3IG model with fitted shift may have been a better choice. The difference models (IG-IG and IG-G) were only a good choice if fit of the tail was most important.

[Figure omitted. See PDF.]

Fig 2. QQ plots for Subject S8 from the awake and at rest cohort.

https://doi.org/10.1371/journal.pcbi.1009099.g002

[Figure omitted. See PDF.]

Fig 3. Rescaled QQ plots for Subject S8 from the awake and at rest cohort.

https://doi.org/10.1371/journal.pcbi.1009099.g003

[Figure omitted. See PDF.]

Fig 4. QQ plots for Subject P10 from the propofol sedation cohort.

https://doi.org/10.1371/journal.pcbi.1009099.g004

[Figure omitted. See PDF.]

Fig 5. Rescaled QQ plots for Subject P10 from the propofol sedation cohort.

https://doi.org/10.1371/journal.pcbi.1009099.g005

Finally, the parameter values (Tables 4–5) indicate that the progressive simplifications are logical, from the difference of two inverse Gaussians to the difference between an inverse Gaussian and a Gaussian, and then to a 3-parameter inverse Gaussian with a location shift. Across all subjects, in the IG-IG model (Model 1), the parameters of the second inverse Gaussian indicate that the ratio of is very high, which approaches a Gaussian distribution [35]. Similarly, in the IG-G model (Model 2), across all subjects, the standard deviation of the Gaussian is very small, approaching a simple point shift in the mean.

[Figure omitted. See PDF.]

Table 4. Fitted parameter values for Models 1–4 for awake and at rest cohort.

https://doi.org/10.1371/journal.pcbi.1009099.t004

[Figure omitted. See PDF.]

Table 5. Fitted parameter values for Models 1–4 for propofol sedation cohort.

https://doi.org/10.1371/journal.pcbi.1009099.t005

Discussion

EDA consists of two simultaneous components or processes at different timescales, the tonic and phasic. Within phasic EDA, there are two sources of relevant physiological characteristics, the timing (temporal information) and size of pulses (amplitudes). In our previous work, we developed a point process model for the temporal information [23,24]. In this paper, we present a model for pulse amplitudes.

We used EDA data from two subject cohorts, a set of healthy volunteers while awake and at rest and another set of healthy volunteers under controlled propofol sedation, to verify our hypothesis that the pulse amplitudes in EDA were characterized by highly regular statistical structure. This statistical structure was consistent with the integrate-and-fire physiology that describes sweat gland function and accounted for the effect of the background, in addition to stimulus intensity, on amplitude. We fitted eight models to the pulse amplitudes in EDA, four of which were simplifications to varying degrees of our hypothesis that each pulse can be modeled as the difference of inverse Gaussians. We quantified the goodness-of-fit with three methods: AIC, and mean and maximum distances from the 45-degree line on the QQ plot. Together, we showed that the model fits were consistent with not only integrate-and-fire sweat gland physiology, but also the combined effects of varying stimulus intensity and EDA background on the dynamics of generated pulses.

The different simplifications of our hypothesis each emphasized fitting different parts of the distribution. For example, the two difference models, Models 1 and 2 (IG-IG and IG-G), both prioritized fitting the right tail of the distribution (the largest pulses) by sacrificing some of the fit near the mode of the distribution. In contrast, the SIG model (Model 4) prioritized fitting the mode of the distribution. The 3IG with fitted shift model balanced both. In the visualization, this can be seen as the varying slope of the QQ-plot relative to the 45-degree line. Going from the most simplified (SIG) to the least simplified (IG-IG) model (from Model 4 to Model 1), each is affected progressively more in terms of slope by the largest pulses. This occurs because the additional parameters, whether the location shift or those of the subtracted distribution, allowed the model to better tailor the fit of the tail.

There were some interesting differences in the performances of the models between the two subject cohorts. In the awake and at rest cohort, the SIG model (Model 4) performed best in terms of AIC and mean distance from the diagonal, while the other simplification models (Models 1–3) all performed well by maximum distance from the diagonal (the 3IG with fitted shift perhaps doing the best). In contrast, in the propofol sedation cohort, the SIG (Model 4) was the best only according to AIC, while the 3IG with fitted shift was the best in terms of both mean and maximum distance from the diagonal. This may reflect a difference in dynamics between both cohorts. Perhaps the more simplified model performed better in the awake and at rest cohort because there were fewer changing dynamics when subjects were largely at rest compared to a changing concentration of drug with known autonomic effects. Or alternatively, perhaps a longer duration of data in the propofol sedation cohort contained more dynamics that required additional complexity in the model.

When examining them together as a framework, the simplification models performed consistently and robustly across both subject cohorts, even though the data were collected under different conditions, using different equipment, and by different researchers. The simplification models were able to capture the structure of pulse amplitude information more successfully than other models across both cohorts, while still being flexible in the degree of complexity of the model as required by the condition under study. The consistent performance of this framework across disparate study conditions and parameters supports the driving physiological model about sweat gland function underlying the framework. If only one subject cohort had been included, there would remain a question of whether the model fits were the result of some specific property of the study condition, equipment, or study design. The inclusion of two different subject cohorts lends support to the physiological validity of the framework.

The result of this study creates a direct link between the physiology of sweat glands and the statistical structure of the pulse amplitude data collected at the skin surface. The most detailed of existing models of EDA are founded in signal processing methods alone, are computationally complex, and make assumptions about pulse amplitudes out of necessity [3–20]. However, looking to the physiology provided a principled framework by which to propose low-order models for pulse amplitudes (maximum of 4 parameters) that account for the effects of both stimulus and background. This result has implications for understanding and tracking the sympathetic component of the autonomic nervous system in a more meaningful way, including both temporal and amplitude information from EDA.

In future work, we will use this result to robustly and accurately capture the valuable physiological characteristics from both the timing and amplitude of pulses in any EDA dataset. We will study the dynamics of both the timing and amplitudes of pulses over time, applying history dependent inverse Gaussian models like those developed for heart rate variability [44–47] and methods for marked point processes [36,37]. We will also study EDA in other contexts, such as during sleep, with pain, and under general anesthesia. Eventually, these methods will have both clinical and non-clinical applications, such as in emotional state and stress detection [18–20]. Our findings provide a principled, physiologically based approach for extending EDA analyses to these more complex and important applications.

Supporting information

S1 Appendix. Additional figures for all subjects.

Table A in S1 Appendix. Model fit results for awake and at rest cohort for Models 5–8. Bold and yellow background indicates best model according to AIC, bold and orange background indicates best model according to mean distance, bold and green background indicates best model according to maximum distance. Table B in S1 Appendix. Model fit results propofol sedation cohort for Models 5–8. Bold and yellow background indicates best model according to AIC, bold and orange background indicates best model according to mean distance, bold and green background indicates best model according to maximum distance. Fig A in S1 Appendix. QQ plots for Subject S1 from the awake and at rest cohort. Fig B in S1 Appendix. Rescaled QQ plots for Subject S1 from the awake and at rest cohort. Fig C in S1 Appendix. QQ plots for Subject S2 from the awake and at rest cohort. Fig D in S1 Appendix. Rescaled QQ plots for Subject S2 from the awake and at rest cohort. Fig E in S1 Appendix. QQ plots for Subject S3 from the awake and at rest cohort. Fig F in S1 Appendix. Rescaled QQ plots for Subject S3 from the awake and at rest cohort. Fig G in S1 Appendix. QQ plots for Subject S4 from the awake and at rest cohort. Fig H in S1 Appendix. Rescaled QQ plots for Subject S4 from the awake and at rest cohort. Fig I in S1 Appendix. QQ plots for Subject S5 from the awake and at rest cohort. Fig J in S1 Appendix. Rescaled QQ plots for Subject S5 from the awake and at rest cohort. Fig K in S1 Appendix. QQ plots for Subject S6 from the awake and at rest cohort. Fig L in S1 Appendix. Rescaled QQ plots for Subject S6 from the awake and at rest cohort. Fig M in S1 Appendix. QQ plots for Subject S7 from the awake and at rest cohort. Fig N in S1 Appendix. Rescaled QQ plots for Subject S7 from the awake and at rest cohort. Fig O in S1 Appendix. QQ plots for Subject S9 from the awake and at rest cohort. Fig P in S1 Appendix. Rescaled QQ plots for Subject S9 from the awake and at rest cohort. Fig Q in S1 Appendix. QQ plots for Subject S10 from the awake and at rest cohort. Fig R in S1 Appendix. Rescaled QQ plots for Subject S10 from the awake and at rest cohort. Fig S in S1 Appendix. QQ plots for Subject S11 from the awake and at rest cohort. Fig T in S1 Appendix. Rescaled QQ plots for Subject S11 from the awake and at rest cohort. Fig U in S1 Appendix. QQ plots for Subject P1 from the propofol sedation cohort. Fig V in S1 Appendix. Rescaled QQ plots for Subject P1 from the propofol sedation cohort. Fig W in S1 Appendix. QQ plots for Subject P2 from the propofol sedation cohort. Fig X in S1 Appendix. Rescaled QQ plots for Subject P2 from the propofol sedation cohort. Fig Y in S1 Appendix. QQ plots for Subject P3 from the propofol sedation cohort. Fig Z in S1 Appendix. Rescaled QQ plots for Subject P3 from the propofol sedation cohort. Fig AA in S1 Appendix. QQ plots for Subject P4 from the propofol sedation cohort. Fig AB in S1 Appendix. Rescaled QQ plots for Subject P4 from the propofol sedation cohort. Fig AC in S1 Appendix. QQ plots for Subject P5 from the propofol sedation cohort. Fig AD in S1 Appendix. Rescaled QQ plots for Subject P5 from the propofol sedation cohort. Fig AE in S1 Appendix. QQ plots for Subject P6 from the propofol sedation cohort. Fig AF in S1 Appendix. Rescaled QQ plots for Subject P6 from the propofol sedation cohort. Fig AG in S1 Appendix. QQ plots for Subject P7 from the propofol sedation cohort. Fig AH in S1 Appendix. Rescaled QQ plots for Subject P7 from the propofol sedation cohort. Fig AI in S1 Appendix. QQ plots for Subject P8 from the propofol sedation cohort. Fig AJ in S1 Appendix. Rescaled QQ plots for Subject P8 from the propofol sedation cohort. Fig AK in S1 Appendix. QQ plots for Subject P9 from the propofol sedation cohort. Fig AL in S1 Appendix. Rescaled QQ plots for Subject P9 from the propofol sedation cohort. Fig AM in S1 Appendix. QQ plots for Subject P11 from the propofol sedation cohort. Fig AN in S1 Appendix. Rescaled QQ plots for Subject P11 from the propofol sedation cohort.

https://doi.org/10.1371/journal.pcbi.1009099.s001

(PDF)

S2 Appendix. Detailed anatomy and physiology.

https://doi.org/10.1371/journal.pcbi.1009099.s002

(PDF)

S3 Appendix. Derivation of method of moments estimates for Models 1 and 2.

https://doi.org/10.1371/journal.pcbi.1009099.s003

(PDF)

Acknowledgments

We would like to thank the MIT Clinical Research Center Staff.

Citation: Subramanian S, Purdon PL, Barbieri R, Brown EN (2021) Elementary integrate-and-fire process underlies pulse amplitudes in Electrodermal activity. PLoS Comput Biol 17(7): e1009099. https://doi.org/10.1371/journal.pcbi.1009099

References

1. Boucsein W. Electrodermal Activity. New York: Springer; c2012.

2. Posada-Quintero HF, Chon KH. Innovations in electrodermal activity data collection and signal processing: a systematic review. Sensors. 2020 Jan;20(479). pmid:31952141

3. Benedek M, Kaernbach C. Decomposition of skin conductance data by means of nonnegative deconvolution. Psychophysiology. 2010;47:647–658. pmid:20230512

4. Benedek M, Kaernbach C. A continuous measure of phasic electrodermal activity. J. Neurosci. Methods. 2010;190:80–91. pmid:20451556

5. Faghih R, Stokes PA, Marin M-F, Zsido RG, Zorowitz S, Rosenbaum BL, et al. Characterization of fear conditioning and fear extinction by analysis of electrodermal activity. Proc. 37th IEEE International Conf on Eng in Biol and Med (EMBC). 2015 Aug.

6. Alexander DM, Trengove C, Johnston P, Cooper T, August JP, Gordon E. Separating individual skin conductance responses in a short interstimulus-interval paradigm. J. Neurosci. Methods. 2005;146:116–123. pmid:15935228

7. Greco A, Valenza G, Lanata A, Scilingo E, Citi L. A convex optimization approach to electrodermal activity processing. IEEE Transactions on Biomed. Eng. 2016;63:797–804. pmid:26336110

8. Hernando-Gallego F, Luengo D, Artes-Rodriguez A. Feature extraction of galvanic skin responses by nonnegative sparse deconvolution. IEEE Journal of Biomedical and Health Informatics. 2018 Sep;22(5):1385–1394. pmid:29990244

9. Amin MR, Faghih RT. Sparse deconvolution of electrodermal activity via continuous-time system identification. IEEE Transactions on Biomedical Engineering. 2019 Sep;66(9):2585–2595. pmid:30629490

10. Jain S, Oswal U, X KS, Eriksson B, Haupt J. A compressed sensing based decomposition of electrodermal activity signals. 2017 Sep;64(9):2142–2151.

11. Bari DS, Aldosky HYY, Tronstad C, Kalvoy H, Martinsen OG. Electrodermal responses to discrete stimuli measured by skin conductance, skin potential, and skin susceptance. Skin Res Technol. 2018 Feb;24(1):108–116. pmid:28776764

12. Posada-Quintero HF, Florian JP, Orjuela-Canon AD, Aljama-Corrales T, Charleston-Villalobos S, Chon KH. Power spectral density analysis of electrodermal activity for sympathetic function assessment. Ann Biomed Eng. 2016 Oct;44(10):3124–3135. pmid:27059225

13. Posada-Quintero HF, Florian JP, Orjuela-Canon AD, Chon KH. Highly sensitive index of sympathetic activity based on time-frequency spectral analysis of electrodermal activity. Innovative Methodology. 2016 Sep;311(3):R582–R591. pmid:27440716

14. Ghiasi S, Grecol A, Nardelli M, Catrambonel V, Barbieri R, Scilingo EP, Valenza G. A new sympathovagal balance index from electrodermal activity and instantaneous vagal dynamics: a preliminary cold pressor study. Annu Int Conf IEEE Eng Med Biol Soc. 2018 Jul;3068–3071. pmid:30441042

15. Posada-Quintero HF, Reljin N, Mills C, Mills I, Florian JP, VanHeest JL, et al. Time-varying analysis of electrodermal activity during exercise. PLoS One. 2018 Jun;13(6):e0198328. pmid:29856815

16. Bari DS, Aldosky HYY, Tronstad C, Kalvoy H, Martinsen OG. Electrodermal activity responses for quantitative assessment of felt pain. J Electr Bioimpedance. 2018 Dec;9(1):52–58. pmid:33584921

17. Bari DS. Psychological correlates of nonspecific electrodermal responses. J Electr Bioimpedance. 2019 Dec;10(1):65–72. pmid:33584885

18. Storm H, Myre K, Rostrup M, Stokland O, Lien MD, Raeder JC. Skin conductance correlates with perioperative stress. Acta Anaesthesiol. Scand. 2002;46:887–895. pmid:12139547

19. Storm H, Shafiei M, Myre K, Raeder J. Palmar skin conductance compared to a developed stress score and to noxious and awakening stimuli on patients in anaesthesia. Acta Anaesthesiol. Scand. 2005;49:798–803. pmid:15954962

20. Sano A, Picard R, Stickgold R. Quantitative analysis of wrist electrodermal activity during sleep. Int. J. Psychophysiology. 2014;94:382–389. pmid:25286449

21. Wickramasuriya DS, Qi C, Faghih RT. A state-space approach for detecting stress from electrodermal activity. 40th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC). 2018 Jul;3562–3567

22. Wickramasuriya DS, Faghih RT. A marked point process filtering approach for tracking sympathetic arousal from skin conductance. IEEE Access. 2020;8:68499–68513.

23. Subramanian S, Barbieri R, Brown EN. Point process temporal structure characterizes electrodermal activity. PNAS. 2020 Oct;117(42):26422–26428. pmid:33008878

24. Subramanian S, Barbieri R, Brown EN. A point process characterization of electrodermal activity. Proc. 40th IEEE International Conf on Eng in Biol and Med (EMBC). 2018 Jul.

25. Subramanian S, Purdon PL, Barbieri R, Brown EN. A model-based approach for pulse selection from electrodermal activity. IEEE Transactions on Biomedical Engineering. 2021. pmid:33822719

26. Subramanian S, Barbieri R, Brown EN. A systematic method for preprocessing and analyzing electrodermal activity. Proc. 41st IEEE International Conf on Eng in Biol and Med (EMBC). 2019 Jul.

27. Kunimoto M, Kirno K, Elam M, Karlsson T, Wallin BG. Non-linearity of skin resistance response to intraneural electrical stimulation of sudomotor nerves. Acta Physiol. Scand. 1992;146:385–392. pmid:1481693

28. Kirno K, Kunimoto M, Lundin S, Elam M, Wallin BG. Can galvanic skin response be used as a quantitative estimate of sympathetic nerve activity in regional anesthesia? Anesth. Analg. 1991;73:138–142. pmid:1854028

29. Wallin BG, Charkoudian N. Sympathetic neural control of integrated cardiovascular function: Insights from measurement of human sympathetic nerve activity. Muscle & Nerve. 2007;36:595–614. pmid:17623856

30. Sonner Z, Wilder E, Heikenfeld J, Kasting G, Beyette F, Swaile D, et al. The microfluidics of the eccrine sweat gland, including biomarker partitioning, transport, and biosensing implications. Biomicrofluidics. 2015;9.

31. Ogawa T, Bullard RW. Characteristics of subthreshold sudomotor neural impulses. Journal of Applied Physiology. 1972;33(3):300–305. pmid:5056204

32. Kunimoto M, Kirno K, Elam M, Wallin BG. Neuroeffector characteristics of sweat glands in the human hand activated by regular neural stimuli. Journal of Physiology. 1991 Oct;442(1):391–411. pmid:1798033

33. Schrodinger E. Zur theorie der fall-und steigversuche an teilchen mit brownscher bewegung. Physikalische Zeitschrift. 1915;16.

34. Gerstein G, Mandelbrot B. Random walk models for the spike activity of a single neuron. Biophys. J. 1964;4:41–68. pmid:14104072

35. Chhikara R, Folks J. The Inverse Gaussian Distribution: Theory, Methodology, and Applications. New York: Marcel Dekker; c1989.

36. Daley D, Vere-Jones D. An Introduction to the Theory of Point Processes. New York: Springer; c2003.

37. Daley D, Vere-Jones D. An Introduction to the Theory of Point Processes: Volume II: General Theory and Structure. New York: Springer; c2007.

38. “Neurofeedback Expert System”, Thought Technology Ltd, https://thoughttechnology.com/neurofeedback-expert-system/, accessed 1/6/21.

39. Purdon PL, Pierce ET, Mukamel EA, Prerau MJ, Walsh JL, Wong KFK, et al. Electroencephalogram Signatures of Loss and Recovery of Consciousness from Propofol. PNAS. 2013;110: E1142–1151. pmid:23487781

40. “BedMasterEx”, Anandic Medical Systems, https://www.bedmaster.net/en/products/bedmasterex, accessed 1/6/21.

41. Subramanian S, Barbieri R, Brown EN. Electrodermal activity of healthy volunteers while awake and at rest (version 1.0). PhysioNet. 2020 Aug. https://doi.org/10.13026/arty-2540

42. Pawitan Y. In All Likelihood. Oxford: Clarendon Press; c2013.

43. Bowman KO, Shenton LR. Estimation: Method of Moments. Encyclopedia of Statistical Sciences. 2006.

44. Barbieri R, Matten EC, Alabi AA, Brown EN. A point-process model of human heartbeat intervals: new definitions of heart rate and heart rate variability. Am. J. Physiol. Heart Circ. Physiol. 2005 Jan;288(1):H424–435. pmid:15374824

45. Barbieri R, Brown E. Analysis of heartbeat dynamics by point process adaptive filtering. IEEE Transactions on Biomed. Eng. 2006;53:4–12. pmid:16402597

46. Barbieri R, Brown E. Application of dynamic point process models to cardiovascular control. Biosystems. 2008;93:120–125. pmid:18515000

47. Valenza G, Akeju O, Pavone KJ, Citi L, Hartnack KE, Sampson A, et al. Instantaneous monitoring of heart beat dynamics during anesthesia and sedation. Journal of Computational Surgery. 2014;1(13).