Point process temporal structure characterizes electrodermal activity

To develop our statistical models of EDA activity we review the anatomy and physiology of sweat production in the skin (1). Each sweat gland consists of three parts: the dermal gland, the duct that connects the gland to the skin surface, and the pore where the duct opens to the skin (Fig. 1A; 1). The dermal portion of the gland is innervated by sudomotor nerves, a part of the peripheral nervous system that is predominantly under sympathetic control (1). Sweat is produced in the gland by sympathetically induced abrupt increases in spiking activity (action potentials) in the sudomotor nerves. These electrical events are called sudomotor bursts. Sweat produced in response to these bursts accumulates in the duct. Once the duct is full, it pushes open the pore and spills onto the skin. The sweat either evaporates or is reabsorbed through the walls of the duct. With the duct now empty, the accumulation process begins again (1, 11, 12).

Sweat on the skin’s surface increases its electrical conductance (inverse of resistance) because its salt content facilitates the propagation of electrical currents. Electrical conductance across the skin can be measured in a standard fashion by placing two electrodes on either the palm or fingers, applying a constant voltage, and measuring the current (1). The pulsatile effects of the sudomotor bursts measured at the skin are termed galvanic skin responses (GSRs). The second-to-second changes in skin sweating measured as second-to-second changes in skin conductance are termed EDA (1).

EDA measurements have two distinct components: tonic activity and phasic activity. Tonic activity represents generally ongoing EDA that reflects the background state of the EDA. The phasic activity reflects primarily the GSR or pulsatile events.

We investigate four statistical models to characterize the times between GSR events recorded during EDA measurements: generalized inverse Gaussian, lognormal, gamma, and exponential. The choice of the generalized inverse Gaussian model follows directly from the physiology described previously. The phasic component of the EDA measurements is comprised of pulsatile events. The times between these pulse events are governed by sympathetic stimulation of the glands, the sweat production and its accumulation in the glands, sweat release on the skin, and sweat reabsorption and evaporation. We represent this four-step sequence of sweat accumulation in the gland and its release onto the skin as an integrate-and-fire process defined by a Gaussian random walk with drift diffusion (Fig. 1B). It is well known that the times between firing events for this elementary diffusion process obey an inverse Gaussian probability model (13–15). Therefore, we chose as the first model the generalized form of the inverse Gaussian probability density defined for an interpulse interval time $x>0$ aswhere $-\infty <\lambda <\infty$, $\psi \ge 0,$ $\chi \ge 0$, and $K_{\lambda }$ is the modified Bessel function of the third kind with index $\lambda$ (16). The generalized inverse Gaussian is a flexible family of probability models that includes the inverse Gaussian model as a special case $\lambda =-\frac{1}{2}$. For $\lambda \le 0$, the generalized inverse Gaussian probability density gives diffusion (integrate-and-fire) models other than the inverse Gaussian. For$\hspace{1em}\lambda >0$, the generalized inverse Gaussian density gives nondiffusion (non–integrate-and-fire) models. As $\chi \to 0$, the generalized inverse Gaussian yields the special case of a gamma probability model with parameters $\lambda >0$ and $2^{-1}\psi \ge 0.$ When $\lambda =1$ as $\chi \to 0$, the generalized inverse Gaussian becomes the exponential probability density (16). Hence, by fitting the generalized inverse Gaussian model we can study systematically how well EDA interpulse interval time series can be characterized by a broad class of integrate-and-fire and non–integrate-and-fire models.

To broaden the class of alternatives to integrate-and-fire models, we also consider the lognormal distribution. Its probability density function is defined aswhere $-\infty <\mu <\infty$ and $\sigma >0.$ The lognormal has been widely used as an empirical model to characterize interevent distributions in point process time series (16).

We postulate that under stable, or approximately stable, background conditions the inverse Gaussian model will provide an accurate description of the inter-GSR (interpulse) event times as the integrate-and-fire model is a plausible description of skin sweat production. This assumes that the group activity of the sweat glands shows a coherent behavior. If there is variation in the behavior of the individual glands, even when the background state is stable, then it is possible that the elementary inverse Gaussian model may not apply.

There may be a mixture of inverse Gaussian models which could produce either longer or shorter interevent intervals relative to a single inverse Gaussian model. These differences could be manifested by either heavier or lighter tails relative to an inverse Gaussian model. To model heavy tail distributions, we consider the lognormal and the gamma distributions, whereas to model the light tail distributions, we consider the gamma and the exponential distributions. The gamma is flexible in that it can have either a heavier or a lighter tail than the inverse Gaussian depending on the values of the parameters (17). Although the exponential distribution is a special case of the gamma, it represents the null model of an underlying Poisson process governing GSR-event production.

We quantify the tail behavior of a probability density $f\left(x\right)$ by evaluating the settling rate, which is defined aswhere $F\left(x\right)$ is the corresponding cumulative distribution function. The settling rate is the limit of the hazard function as $x$ tends to infinity (17). By this definition, a distribution is commonly classified as light, medium, or heavy tailed based on whether the settling rate is infinite, positive but finite, or zero, respectively. The more slowly (rapidly) the tail of a distribution settles, the heavier (lighter) the tail. We can divide our set of probability models into four medium-tailed (exponential, gamma, inverse Gaussian, and generalized inverse Gaussian) and one heavy-tailed distribution (lognormal). The four medium-tailed distributions are distinguished relative to each other based on their parameter values, which determine the respective settling rates (17, 18).

Our experimental protocol was approved by the Massachusetts Institute of Technology Institutional Review Board, and all subjects provided written informed consent. We collected EDA data from 12 healthy volunteers (six men) between the ages of 22 and 34 while awake and at rest (19). Electrodes were connected to the second most distal phalange of the second and fourth digits of each subject’s nondominant hand. Approximately 1 h of EDA data were collected at 256 Hz. 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 use the instrumented hand. We assumed skin and ambient temperature were constant for the duration of the experiment. One subject’s data were not included in the analysis because we learned, after the data collection, that this individual occasionally experienced a Raynaud’s type phenomenon. This would affect the quality of the EDA data. Data from the remaining 11 subjects were analyzed using MATLAB 2017a.

Preprocessing consisted of two steps: 1) detecting and removing artifacts and 2) isolating the phasic component. Artifact detection was done based on the derivative of the time series since large rapid changes are physiologically impossible for skin conductance. Artifact removal was done in two parts, first correcting for artifact-related large magnitude changes in the remainder of the signal and then interpolating the few seconds around the artifact itself. Then a low-pass FIR filter was used to estimate and remove the slow-moving tonic component of the signals, thereby, isolating the phasic component. The data preprocessing is described in further detail in Subramanian et al. (20).

Since there are underlying tonic fluctuations, absolute pulse amplitude alone is insufficient to extract pulses reliably. Therefore, we computed locally adjusted amplitudes for all detected peaks using the MATLAB function findpeaks. The findpeaks algorithm computes a prominence or relative amplitude for each peak, by adjusting the amplitude of each peak 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. With this method, a peak with small absolute amplitude can be rewarded in its prominence value if it is in a region of data with low activity. We then used a threshold on this prominence value to account for the varying baseline, instead of on the absolute amplitude.

We used a prominence threshold of 0.005 to extract peaks across all subjects, unless this resulted in too few or too many pulses for each hour-long recording. This was primarily verified by visual inspection of extracted pulses, as well as rough estimates of 60 and 360 pulses as generous bounds for what should be expected of 1 h of data at rest with little external stimulation. This corresponds to one pulse every 10 to 60 s on average (4, 21–26). In the case that too few pulses were extracted, we gradually reduced the prominence threshold by 0.001 until the number of pulses exceeded 100 and no obvious pulses were missed by visual inspection (verified by three different viewers independently). In the case too many pulses were extracted, we gradually increased the prominence threshold by 0.001 until the number of pulses was fewer than 360 and the extracted pulses were distinguishable from sensor noise by visual inspection (verified by three different viewers independently). If changing the threshold in increments of 0.001 resulted in drastic changes in the number of pulses each time, we reduced the increment to 0.0005. Although not fully automated at this stage, the design of this method to extract pulses attempted to take into account the wide variation in baseline levels of EDA activity seen across subjects. The extracted peaks included smaller peaks that other methods would generally ignore as noise. However, we chose to include them in the analysis.

First, we restricted the parameter space to $\lambda \le 0$ for the generalized inverse Gaussian model to find the best integrate-and-fire model for each data series. To allow for statistical improvements by capturing deviations from the integrate-and-fire models, we fit four other probabilities models to each data series: generalized inverse Gaussian non–integrate-and-fire $\left(\lambda >0\right)$, lognormal, gamma, and exponential models. We fit all models by maximum likelihood (27, 28). We assessed goodness-of-fit by using KS plots and by computing Akaike’s information criterion (AIC), defined aswhere $f\left({\widehat{\theta }}_{ML}\right)$ is the likelihood evaluated at the maximum likelihood parameter estimates and $p$ is the number of parameters. A lower AIC indicates a better fit.

A KS plot compares the rescaled quantiles from the fit of the estimated probability model with the quantiles of an exponential distribution with rate 1 using the time-rescaling theorem (29). This theorem states that any point process can be rescaled to an exponential distribution with rate 1 using its conditional intensity (hazard) function. The KS distance computes the maximum distance between the quantiles of the rescaled data and the uniform distribution, which is a simple transform of an exponential distribution with rate 1. A smaller KS distance indicates that the model is more similar to the structure observed in the data. We computed 95% confidence intervals (5% significance cutoffs) for the KS plot and compared the KS distances across models (30). A KS distance that is within (outside) the 95% confidence intervals suggests that the model offers (fails to offer) a reasonably accurate description of the data.

We also compared the models using a tail behavior analysis, in which the settling rates of the models were compared to determine the heaviness of the tails of the distributions. We hypothesized that each EDA interpulse interval distribution results from the activity of multiple sweat glands. This could lead to deviations from the inverse Gaussian with slightly heavier or lighter tails. These deviations can be captured statistically by right-skewed models such as the lognormal, generalized inverse Gaussian and gamma that together allow for more flexibility in tail behavior.

Fig. 1C shows an example of an excerpt of extracted pulses for Subject S07. This includes pulses large enough to be used in most analyses as well as those that are much smaller and usually either smoothed out or ignored as noise. We included both types of pulses for all subjects and did not distinguish between them. The majority of subjects showed appreciable fluctuations in the tonic component across time (SI Appendix, Fig. S1). Across the 11 subjects, the total number of pulses in the 1-h time window ranged between 97 and 348, including the distantly spaced smaller pulses (SI Appendix, Fig. S2). The final prominence thresholds used ranged between 0.0025 and 0.023.

For all 11 subjects, the optimal generalized inverse Gaussian diffusion model was an inverse Gaussian, indicated by λ=−0.5 in Table 1. This inverse Gaussian was always within the significance cutoff according to KS distance, indicating that it provides an accurate description of the data and supporting our hypothesis that the pulsatile sweat release events of EDA can be modeled as an integrate-and-fire process. When allowing for deviations from the inverse Gaussian, for 9 of the 11 subjects, one or more of the lognormal, gamma, or generalized inverse Gaussian nondiffusion models was able to improve statistically on the fit of the inverse Gaussian (Tables 2 and 3; Fig. 2; and SI Appendix, Figs. S3–S7). The exponential was never a better fit than the inverse Gaussian, and was only within the significance cutoff for 4 of the 11 subjects. This suggests that for the majority of subjects, the exponential model did not offer an accurate description of the data.

Both AIC and KS distance were in agreement for the best fit models for the data of 8 of the 11 subjects. This corresponded to 1 generalized inverse Gaussian diffusion model (S08), 1 generalized inverse Gaussian nondiffusion model (S10), 4 lognormal (S02, S03, S09, and S11), and 2 gamma (S04 and S05) models. For S01 and S06, the AICs and KS distances suggested different best fit distributions between inverse Gaussian and lognormal. For S07, they identified different best fit distributions between lognormal and gamma. It is reasonable to expect that the results from AIC and from KS distance will not match exactly since different metrics are intended to capture different aspects of model fits. However, the fact that they agree across the majority of subjects reinforces that there is specific statistical structure in the data that can be captured with a parsimonious model.

The second phase of comparing the models was analyzing the tail behavior (Table 4) using the settling rates estimated for each of the five models for all subjects. Across all 11 subjects, the lognormal always had the heaviest tail since it is commonly classified as a heavy-tailed distribution. The generalized inverse Gaussian integrate-and-fire model had the next heaviest tail, indicated by the next smallest settling rate. The lognormal was the only other class of models besides the generalized inverse Gaussian diffusion models that was within the significance cutoff for all subjects. This suggests that deviations from the inverse Gaussian tended toward heavier tails that were best captured by the lognormal.

Among the remaining models, the generalized inverse Gaussian non–integrate-and-fire models always had lighter tails than the integrate-and-fire models. The tails of the gamma were even lighter. Omitting the exponential, since it fit poorly in the majority of cases, the remaining models—lognormal, generalized inverse Gaussian (diffusion and nondiffusion), and gamma—together provide a systematic framework to: 1) evaluate the presence of inverse Gaussian structure in EDA using diffusion models; and 2) enhance statistical descriptions using models capable of capturing a range of tail behavior properties.

Among the three subjects for whom the AIC and KS distance disagreed on the best fit model (subjects S01, S06, and S07), there was only one case (S07) in which this disagreement predicted very different tail behavior. In this case, the AIC and KS distance were divided between the lognormal and the gamma as the best model, one of which predicted a very heavy tail while the other predicted a very light tail. However, the tail of the interpulse interval distribution may have been relatively underrepresented for this subject, so the fit was primarily based on the shape near the mode, which both the lognormal and gamma distributions fit well. With a data collection period longer than 1 h, we would most certainly arrive at a more accurate description of the tail behavior.