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

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.
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][4][5][6][7][8][9][10][11][16][17][18][19][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][28][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.

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][4][5][6][7][8][9][10][11][16][17][18][19][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][28][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][28][29][30][31][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 integrateand-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][34][35].
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), 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-andfire 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), 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 as where f ðŷ ML Þ 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.
We also plotted QQ plots and rescaled QQ plots. For rescaled QQ plots, both model and empirical quantile values were rescaled by q rescaled ¼ 1 À expðÀ qÞ: Rescaled 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).

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 Exponential f x; l ð Þ ¼ le À lx https://doi.org/10.1371/journal.pcbi.1009099.t001

PLOS COMPUTATIONAL BIOLOGY
Elementary integrate-and-fire process underlies pulse amplitudes in Electrodermal activity 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). 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 Table 2. Model fit results for awake and at rest cohort for Models 1-4. 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.

AIC
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 l 2 m 2 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.

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

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

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  Table 4. Fitted parameter values for Models 1-4 for awake and at rest cohort. 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.

SIG (Model 4) 3IG with fitted shift (Model 3) IG-G (Model 2) IG-IG (Model 1)
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][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][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][45][46][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][19][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