Robust Inference of Autonomic Nervous System Activation Using Skin Conductance Measurements: A Multi-Channel Sparse System Identification Approach

The autonomic nervous system (ANS) stimulates various sweat glands for maintaining body temperature as well as in response to various psychological events. Variations in skin conductance (SC) measurements due to salty sweat secretion can be used to infer the underlying ANS activity. Recovering both ANS activity and the underlying system from noisy single-channel recordings is challenging. As the same ANS activity drives all the sweat glands throughout the skin, the same information is encoded in different SC recordings. We perform system identification and develop a physiological model for multi-channel SC recordings relating them to ANS activation events. Using a multi-rate formulation, we estimate the number, timings, and amplitudes of ANS activity and the unknown model parameters from multi-channel SC data. We incorporate a generalized-cross-validation-based sparse recovery approach to balance between the sparsity level of the inferred ANS activity and the goodness of fit to the multi-channel SC data. We successfully deconvolve multi-channel experimental auditory stimulation SC data from human participants. We analyze experimental and simulated data to validate the performance of our concurrent deconvolution algorithm; we illustrate that we can recover the ANS activity due to the underlying auditory stimuli. Furthermore, we estimate stress using inferred ANS activity based on multi-channel deconvolution of SC data collected during different driving conditions and at rest. We propose a model for multi-channel SC recordings. Moreover, we develop a multi-channel deconvolution approach to perform robust sparse inference in the presence of noise. The proposed approach could potentially improve stress state estimation using wearables.


I. INTRODUCTION
Electrodermal activity (EDA) refers to any alterations in the electrical characteristics of the skin caused by salty sweat secretion.Hypothalamic control of sweating is primarily intended for thermoregulation of the human body.Apart from thermoregulation, sweating can also take place due to other physiological events including emotional arousal [1], [2].Moreover, variations in skin conductance (SC), a mea-The associate editor coordinating the review of this manuscript and approving it for publication was Kezhi Li .sure of EDA, are highly correlated with emotions and can be used for interpreting emotional dysregulation and disturbances [3], [4].Alterations in SC throughout the different skin regions are regulated by the autonomic nervous system (ANS) [5].The analysis of SC time series assumed to be modulated by the ANS can potentially be used to track the mental health of an individual in order to prevent mental stress-related problems [6].
Walker et al. [7] reported that a large portion of the deaths worldwide are attributable to mental health-related disorders.Regular tracking of problematic patterns of emotional regulation could potentially help prevent psychiatric disorders [8].Physiological signals including electroencephalogram, electrocardiogram, respiration, functional near-infrared spectroscopy [9], and EDA, could be investigated to identify abnormal patterns of emotional regulation [10].Day-to-day monitoring and tracking of emotional regulation require wearable implementations.However, a reliable, noise-robust monitoring system is challenging and yet to be developed.A noise-robust personalized wearable mental health monitoring system could eventually lead to effective monitoring of mental health-related problems [11].
In another context, several recent studies have shown correlations between abnormal SC recordings due to diabetic neuropathy and other different diabetic diseases, in different races [12]- [14].Diabetic neuropathy is a type of nerve injury caused by diabetes [15].Nerves in the legs, feet, and hands are more prone to this type of damages [15].These nerves include the sudomotor nerves that are responsible for delivering the ANS stimulation to the sweat glands.Diabetic neuropathy may cause abnormal SC in different regions of the body due to nerve damage.Appropriate systematic analysis of SC recording from different regions may lead to early diagnosis and prevention of diabetes-related complications.
As the SC measurement can be thought of as the convolution of ANS activity with physiological smoothing kernels that consist of sweat glands, sweat secretion, and evaporation dynamics [16], identification of ANS activity is a deconvolution problem.Many research works have been carried out with different methods for the deconvolution of SC recordings to recover the timings and amplitudes of stimulation and to estimate the underlying physiological parameters with the goal of uncovering emotional states using single channel SC data.Benedek and Kaernbach [17], [18] presented a scheme where a non-negative deconvolution method is utilized to separate single SC time series into discrete compact responses.They have also analyzed SC responses to assess the deviations from the standard SC response shape.Nevertheless, their scheme does not consider sparsity constraint to prevent capturing noise as SC responses.Moreover, they do not consider the individual differences in the modeling of the rise and decay times.They perform their deconvolution for multiple predefined sets of parameters and choose the one that provides the most reasonable fit.
Greco et al. [19] proposed a convex optimization formulation to decompose SC time series into tonic and phasic components.Unlike in [18], they considered a the sparsity condition in neural stimuli from the ANS.They use a fixed regularization parameter for imposing the sparsity constraint.However, finding an accurate sparse solution that can handle inter-and intra-subject variability is challenging with a fixed regularization parameter for automatic processing.In a similar work, Hernando-Gallego et al. [20] proposed a faster decomposition approach to obtain a sparser solution; however, this approach leads to overly sparse solutions compared to the underlying neural stimuli.Faghih [21], Faghih et al. [22], Faghih [23], and Faghih et al. [24] proposed a two-step coordinate descent deconvolution approach to account for individual differences in the physiological system parameters for SC signals and similar physiological system.These approaches have been successfully utilized to characterize hormonal dysfunction in fibromyalgia patients [25] as well as potential neurobiologic substrate for chronic insomnia [26].However, inference of ANS stimulation along with physiological system parameters using single channel SC recording is challenging in the presence of noise.
A substantial amount of research has been carried out with wearable EDA sensors that can provide better insight into how affect and stress interact with daily life [27]- [29].However, wearable sensors often suffer from poor signal quality as well as motion artifacts.Furthermore, developing algorithms for analyzing EDA data in the presence of artifacts and noise is yet to be undertaken.Fusing multi-channel recording with wearable devices to account for the poor signal quality could potentially improve the deconvolution performance.Faghih et al. [30] proposed an approach to include multi-channel hormone time series using a combined statespace model, and then, developed a concurrent deconvolution scheme.
Inspired by the work by Faghih et al. [30] for concurrent deconvolution of the multi-channel hormone time series, we utilize concurrently collected multi-channel SC measurements to robustly infer ANS stimulation.We hypothesize that the changes in SC in different regions of skin are due to the same ANS stimulation.We implement system identification and propose a state-space model that includes the SC recordings from multiple skin regions.Furthermore, we introduce a multi-channel concurrent deconvolution algorithm and analyze SC data using auditory stimulation experimental data.Finally, we use a Bayesian approach [31]- [33] to track a cognitive stress level based on the concurrent deconvolution of SC data collected during different stressful driving conditions.

A. DATASET DESCRIPTION 1) DATASET 1 -AUDITORY STIMULATION
We analyzed the SC responses to loud sounds, simultaneously recorded from palm, fingers and foot data [34].The experiment was conducted for modeling event-related SC responses.The dataset contains SC data recorded from 26 healthy participants from three different skin locations: the hypothenar of the non-dominant hand, the middle phalanx of second and third finger of the dominant hand, and the medial plantar surface of the non-dominant foot.Each participant was provided 20 auditory stimuli.Each stimulus is single white noise burst of 1s duration.Participants were asked to press a foot operated pedal upon hearing the stimuli.A detailed description of the auditory stimulation experiment is given in [35].We discard the data contaminated with heavy motion artifacts and the data having very small SC responses from our study.

2) DATASET 2 -DRIVER STRESS
To assess the performance in tracking stress using deconvolution result, we also analyzed the stress recognition in the automobile drivers dataset [36].This dataset includes recordings for 17 separate driving sessions on a predefined route with highways, toll roads, and city driving.The driving sessions were all conducted during the mid-morning or mid-afternoon with light traffic.Rest periods were included during which the drivers sat in the car with their eyes shut.Since annotations of each portion are not publicly available for all the subjects, we could only use one recording whose approximate timings had to be matched with a figure in [36].We apply the stress tracking algorithm by Wickramasuriya et al. [31] on the deconvolution result from the current work and compare the result with the heuristic approach in [31].

B. MODEL FORMULATION FOR PHASIC SC DECONVOLUTION
The variations in SC, caused by salty sweat secretion, are regulated by sudomotor nerve activity (SMNA) in the ANS.SC data can be interpreted as the summation of a fastvarying component and a slow-varying component.The slowvarying component, known as the tonic component, is mainly intended for thermoregulation of the body.The comparatively fast-varying component is termed as the phasic component.The phasic component represents the activity of ANS, which is a reflection of emotional events.The SC signal can be represented as the summation of tonic and phasic components as follows: where y sc n (t), p n (t), and s n (t) represent the SC signal and its phasic and tonic components for the nth channel, respectively.In the pre-processing stage, we extract the phasic components p n (t) from the SC data y sc n (t) using cvxEDA [19].In this case, we use the default value of the regularization parameter for l 1 -norm minimization in cvxEDA [19].We assume that the extracted phasic components have some Gaussian noise added to them.We consider this as measurement noise.We hypothesize that the same ANS stimulation is responsible for modulating the phasic component of SC data related to different eccrine sweat glands in different regions of skin.We propose a state-space model with n-channel phasic SC data and the single ANS input.Figure 1 shows the complete block diagram of the proposed system model.We model p n (t) as a scaled version of a representative internal state ζ n (t), which refers to the average amount of sweat in the epidermis associated with the nearest sweat gland.We introduce a scaling factor α n to account for the number of sweat glands present per unit area where the nth sensor is placed.Relation between p n (t) and ζ n (t) is as follows, We describe the system dynamics using the following set of differential equations each denoting the kinetics of sweat secretion and evaporation process in sweat glands [16], [24], [37], where ζ n (t) is the internal state which is reflected into y n (t) from the nth channel, ∀n ∈ {1, 2, • • • , χ }.The term β n refers to the delay in the stimuli input for the nth channel.The terms τ r and τ d correspond to the rise and decay times of the SC responses, respectively; these parameters are fixed for all channels.β n can be calculated by taking the crosscorrelation of the y sc n and y sc 1 before deconvolution.The location of the maximum value of the cross-correlation is used to calculate the time lag β n .We assume β 1 = 0 and α 1 = 1 for the reference channel.Under the assumption that the ANS stimuli are sparse as in [24], we represent the input as u(t) = N i=1 q i δ(t − i ) where q i is the amplitude of the impulse in the neural stimuli at time i .
Let x 2n−1 (t) and x 2n (t) = ζ n (t + β n ) be the internal states for ∀n.Similar to [24], [37], the differential equation in (2) for the nth channel can be re-written in state-space form as follows, The corresponding continuous observation equation can be written as follows, where y n (t) is the continuous observation variable and ν n (t) refers to the noise process.In matrix form the state-space model can be written as follows, We can write the equations in ( 6)-( 7) in a state-space form, for all channels, as follows, where and Discretization: Let T u be the sampling frequency of the neural stimuli and T y be the sampling frequency of the phasic SC data for each channel.The timings of the neural impulses can be written as i = iT u ; q i is zero if there is no impulse at the ith instance.Let y n, k be the observed phasic SC for the nth channel at time instance t k = kT y .We can write where ν n, k is the noise associated with the nth channel; ν n, k is modelled as a zero-mean Gaussian random variable.We derive the discrete equivalent of the system, assuming that the input and the states are constant over T u .The discrete version of the neural stimuli can be written as a vector u = [q 1 q 2 • • • q N ] that represents the entire neural stimuli over the duration of SC data.Let = e A c T u , and = T u 0 e A c (T u −ρ) B c dρ to write the discrete state-space form of ( 8)-( 9) as: As neural stimuli and SC measurement have different sampling frequencies, i.e.T y = LT u where L is a positive integer, we let ; the multi-rate system can be represented as follows: where A d and B d are functions of τ = τ r τ d , α, T u , and T y .Let θ = τ α .As the system is causal, we use ( 13)-( 14) to obtain the observation equation for the kth sample: where For the initial condition, we can let similar to the work in [24]. where Therefore, we can write the solution for the observation equation for all the sampled data as follows, Equivalently, we can separately represent the solution for each channel as follows, . . .

NOISE VARIANCE ESTIMATION
We assume that the noise term ν n,k is Gaussian random variable with zero-mean and σ 2 n variance.Therefore, the energy of the noise is distributed over the whole spectral range of 0 to half of the sampling frequency.To obtain noise variances ∀n, we filter the phasic components of the experimental SC signals with a 0.5 Hz high pass filter to remove all the signal components assuming the signals are band-limited to 0.5 Hz.Then, we calculate the variances of the filtered signals to obtain the estimates of the noise variances in the high-frequency region.We interpolate these estimated variances for the whole spectral bandwidth considering the lowfrequency region.This enables us to obtain an assessment of the noise variances σ 2 n , ∀n.

C. MODEL FORMULATION FOR STRESS TRACKING
We assume that the ANS stimulation u for SC responses in SC data is dependent on an internal stress state in the brain.Wickramasuriya et al. [31] modeled the evolution of the stress state as a random walk.SC data is modulated by stress and thus the probability p j of occurrence of a neural stimuli is dependent on the stress state w j .Given the state process w j , the observation model defines the probability of observing a neural stimulus.Let the observation variable s j be a binary variable where s j = 0 denotes no neural stimulus and s j = 1 denotes a single neural stimulus observed at jth time instance.The probability distribution of s j can be described using the Bernoulli distribution.State-space model is as follows, w j = w j−1 + j (state equation) ( 16) P(s j |q j ) = p s j j (1 − p j ) 1−s j (observation equation) (18) where p j is defined by a logistic equation.The observation model in (18) is defined using the Bernoulli distribution.The parameter η represents the probability of observing a neural stimulus by random chance in a bin at the start of the experiment.We calculate η similarly as in [31].We use a bin size of 1 second in u to generate S 1:J .Here J denotes the total number of bins.As the signal sample frequency is also 1 second, number of observation is equal to the number of samples per channel, i.e., J = M .With the observation S 1:J = {s 1 , s 2 , . . ., s J } indicating the presence or absence of neural stimuli, we estimate W = {w 0 , w 1 , w 2 , . . ., w J } and σ 2 for in turn estimating p j , ∀j using the Expectation Maximization (EM) algorithm in [31].A brief description of the EM algorithm is given in Subsection II-D3.

D. ESTIMATION 1) PREPROCESSING
We discard all the participants having very heavy motion artifacts in their data for our analysis.The original recorded signal has a sampling frequency of 100 Hz for dataset 1 and 15.5 Hz for dataset 2. We perform all the preprocessing in the original sampling frequency.In the beginning preprocessing step, we find the large discontinuities in the recorded SC data.We take the difference and the negative of the difference of the raw SC data and detect peaks having a prominence of 0.1.We assume these discontinuities are due to the artifacts, and we choose the prominence parameter to detect the artifacts that are much larger than the phasic responses.For a detected discontinuity, we discard a small patch of 0.5 seconds, while keeping the discontinuity in the middle of the patch.We interpolate the region of the discarded patch in the signal with a spline curve.Then, we perform lowpass filtering on the signal with a 64 order FIR lowpass filter of 3 Hz cut off frequency.As the FIR lowpass filter has linear phase, we also corrected the group delay generated by the FIR filter.Afterward, we apply cvxEDA [19] to separate the tonic component and the phasic component, taking time between knots of the tonic spline function as 6 seconds.All other parameters are kept unchanged from their default value.Next, we calculate the delay parameter β taking cross-correlation between the two phasic components from the first channel and the nth channel for different lags.We take the lag for foot phasic SC data with the maximum correlation.Finally, we resample the signal to 1 Hz sampling frequency.

2) DECONVOLUTION
The sampling period for phasic SC signal and neural stimuli are T y = 1 seconds and T u = 0.25 seconds, respectively.To estimate the system parameters and the neural stimuli, we use solutions in (15), and formulate the following optimization problem imposing the sparsity constraint on u: where χ+1) , and b ∈ R 2(χ+1) .The l p -norm is an approximation of the l 0 -norm (0 < p ≤ 2) [22].The l p -norm regularization parameter λ is chosen to maintain a balance between filtering out the noise and the sparsity level of u [22], [30], [37].We solve the inverse problem of finding a non-negative u with a specific sparsity level using the Focal Underdetermined System Solver (FOCUSS+) algorithm [38].Afterwards, we use the generalized cross-validation (GCV) technique to calculate an appropriate value of the l p -norm regularization parameter λ adaptively similar to the approaches in [22], [24], [30].

G(λ)
5) Iterate until convergence We solve the problem in (19) using the following algorithm: We run the algorithm with several random initialization of system parameters τ r , τ d and α n ∀n.Finally, we choose the estimate that has the minimum value for the cost function in (19).

3) STRESS ESTIMATION
We follow the Expectation Maximization (EM) approach in [31] to estimate the stress state w j .In the E-step of their proposed algorithm, a nonlinear recursive point process forward filter followed by a backward smoother has been implemented to estimate the stress state of the subject.The forward filter estimates the stress state w j|j in the jth bin, given S 1:j , i. e., the observations up to the jth bin.The backward smoother estimates the stress state w j|J in the jth bin, given S 1:J , i. e., all the observations (up to the J th bin).p j|j denotes the probability of a neural stimuli impulse occurring within the jth bin given the observation S 1:j , and p j|J the probability of a neural stimuli impulse at jth bin given all of the data S 1:J .In this case, the value of parameter σ 2 has been determined by maximizing the complete data log-likelihood likelihood estimate at previous iteration.= û(i) and λ(i iii.Set λ = λ(i) (m−1) and θ = θ(i−1) ; solve for û(i) (m) by initializing the optimization problem in (19) .

b: MAXIMIZATION STEP
At the maximization step, the expected value of the complete data log-likelihood is used to select the model parameters for the next iteration as follows [40]: The EM algorithm repeatedly iterates between the E-step and the M-step until convergence.The estimated state w j at each time instance is assumed to be Gaussian distributed w j ∼ N (w j|J , σ j|J ) and we define high arousal index (HAI) as Pr(w j > w T ) similar to [32].The threshold w T is set to subject's mean stress state value across the whole experiment.

A. DATASET 1
We use the proposed algorithm and concurrently deconvolve SC measurements from the middle phalanx of hand and the medial plantar surface of foot collected during an auditory stimulation experiment and recover the underlying stimuli u(t), the corresponding rise time (τ r ) and decay times (τ d ) of SC responses, and the attenuation (α 2 ) at the medial plantar surface of foot with respect to the middle phalanx of hand.Results in Figure 2 show that the proposed algorithm successfully recovers the timings and amplitudes of neural stimuli and the underlying system parameters, i.e. the rise and decay times for two female participants and two male participants.The figures for the deconvolution result for all 12 participants are given in the Appendix.We considered the signal segment from 200 seconds to 400 seconds for our analysis.The multiple correlation coefficient (R 2 ) has been calculated for all twelve reconstructed signals.The high values of R 2 (found to be greater than 0.95 except for male participant 3 and the explanation is given in Section IV) for hand phasic SC data suggest that our proposed algorithm can successfully recover the underlying physiologically plausible ANS stimulation.For foot data, the R 2 values are greater than 0.80 except for the male participant 4 and it also explained in Section IV.In general, the R 2 values from reconstructed foot data are lower compared to the R 2 values from reconstructed hand data.Lower R 2 values for foot data suggest that foot signals are noisier than the hand signals.Table 1 shows all the estimated parameters as well as the R 2 values for all 12 participants.R 2 h and R 2 f in Table 1 correspond to the multiple correlation coefficient (R 2 ) for the hand and the foot channels, respectively.The quantile-quantile plot of the model residuals after deconvolution follows a straight line, suggesting that the residuals are Gaussian distributed.All the quantile-quantile plots for 12 are provided in the Appendix.Figure 3 shows an example of three channel deconvolution.Here we include the third recording from the thenar/hypothenar of the hand in the concurrent deconvolution scheme.We are able to successfully deconvolve the three channel SC data and obtain all the unknown parameters.The R 2 obtained for all three channels are 0.985, 0.937 and 0.965.
To validate our approach, we simulate data using the results from the deconvolution and the parameters obtained are given in Table 1.To simulate noisy data, we added zero-mean Gaussian noise with 30 dB and 20 dB signal to noise ratio (SNR) to the reconstructed phasic SC data corresponding to both hand and foot channels, respectively.These SNR values are chosen to obtain comparable noise levels as in experimental data.Thereafter, we again deconvolve the noisy simulated data to compare the results with the ground truth used to simulate the data.Figure 4 shows results from four simulated data.The deconvolution results for all twelve participants are given in the Appendix.Table 2 shows the estimated parameters, the corresponding estimated errors, and R 2 values.
To further validate our approach, we simulated noisy data using a synthetic input and the model parameters τ r = 0.75, τ d = 4, and α 2 = 0.3.In this case, we have the ground truth for comparison.To simulate noisy data, we added zeromean Gaussian noise with 20 dB and 15 dB SNR to the hand and foot phasic SC data, respectively.These two levels of noise are chosen bacause of the higher levels of noise in foot data.Figure 5 shows the simulated data for both channels.We perform deconvolution on the simulated noisy data.Figure 5 also shows the deconvolution performance of other existing algorithms for comparison, as well as our single channel deconvolution approach.Results using simulated data show that our concurrent deconvolution scheme outperforms existing methods.The last three panels in Figure 5 and the corresponding estimation errors in Table 3 also show how our concurrent deconvolution scheme performs better than our single channel deconvolution approach.Moreover, Figure 6 shows how noise effects the estimation accuracy for τ r and τ d .

B. DATASET 2
Figure 7 shows the stress estimates for the estimated neural stimuli using the heuristic peak detection [31], our single channel deconvolution [32], and our concurrent deconvolution.Figure 7a shows the stress estimation result using the approach proposed in [31] for the driver using heuristic peak detection scheme for detecting neural stimuli u.We perform both multi-channel and single channel deconvolution on a small segment from 1500 seconds to 1700 seconds of the recorded signal.After deconvolving phasic SC segment, we solve the inverse problem of estimating u for entire phasic SC signal using FOCUSS+ and GCV-FOCUSS+ with the θ estimated from the deconvolution step.Then, we bin the estimated u with a bin size of 1 seconds.We use the same stress estimation approach in [31] with the result obtained using our concurrent deconvolution algorithm.

IV. DISCUSSION
Figure 1 shows the overall system block diagram, where the system parameters for different regions are assumed to be same and the attenuation terms 0.01 ≤ α 2 , α 3 , • • • , α χ ≤ 100.However, we have also considered other configurations for the system and tried to formulate optimization problems.We investigate different possibilities with two channel recordings.We have considered the configuration by taking τ the same for both hand (middle phalanx) and foot (medial plantar surface) regions.In this case, the attenuation α 2 = TABLE 1.The estimated model parameters and the squares of the multiple correlation coefficients (R 2 ) for the fits of the experimental skin conductance time series.

FIGURE 4. Estimated Deconvolution of the Simulated Phasic SC Signal:
The panels show the deconvolution result on the simulated data for two female and two male participants, respectively.In each panel, i) the top sub-panel shows the simulated (blue stars) and the estimated (red curve) phasic component corresponding to the middle phalanx of hand; ii) middle panels shows the simulated (blue stars) and estimated (red curve) phasic component corresponding to the the medial planar surface of foot; and iii) bottom sub-panel shows the timings of the simulated ANS activation timings and amplitudes (gray line) and the estimated ANS activation timings and amplitudes (red dashed line).
1 which corresponds to the recording from foot.The R 2 squared fit for the foot is very low and a clear attenuation factor other than α 2 = 1 is observed.Then, we formulate the problem with the configuration, which has two different sets of system response parameters (rise time and decay) τ for both channels while keeping α 2 = 1.This time, the solutions for τ for both hand and foot stagnated at extreme bounds and fit after reconstructions were not good for at least one channel.Next, we have analyzed the data taking τ different and 0.01 ≤ α 2 ≤ 100.In this case, the fits TABLE 2. The estimated model parameters and the squares of the multiple correlation coefficients (R 2 ) for the fits of the simulated skin conductance time series.

FIGURE 5. Performance Comparison of Proposed Concurrent
Deconvolution Approach with Existing Approaches: The panels (i) and (ii) respectively depict the synthetic simulated data with 20 dB and 15 dB noise for hand and foot.The panels (iii) and (iv) show results with LedaLab [18].The panels (v) and (vi) show the results with cvxEDA [19].The panels (vii) and (viii) show the recovered neural stimuli with our single channel deconvolution.The panel (vii) shows the recovered results with our concurrent deconvolution with simulated hand and foot SC data.Gray vertical lines correspond to the ground truth, and red, blue and green vertical lines corresponds to recovered neural stimuli with hand data, foot data, and concurrent deconvolution.
after reconstruction were good for most of the channels but all the solutions for τ for both hand foot regions stagnated in the upper/lower bounds.This is because the optimization formulation has too many degrees of freedom.Therefore, we remove some degrees of freedom by taking τ same for both hand and foot regions.In addition to that, we introduced a delay in the input in the foot sweat glands causing delayed phasic response and successfully, concurrently deconvolved phasic SC signals from hand and foot, and recovered neural stimuli and underlying physiological system parameters.
It is somewhat counter-intuitive that the skin conductance system responses (rise times and decay times) in different regions are similar although there can be dissimilarities in the sweat glands.This is because different rise times and decay times for different skin locations will lead to too many degrees of freedom.It is possible to have many solutions to the system for a given set of observed data.In this work, we assumed the system responses does not have much variance in different skin regions.Moreover, the variability in different skin regions is captured by the tonic component.We remove the tonic component before proceeding for deconvolution steps.In one of our works [41], we have developed a deconvolution scheme to separate the tonic and phasic component iteratively in a coordinate descent manner.It is possible to develop a similar approach for multi-channel case.In this way, the variability of the phasic SC response in different skin locations can be captured.It might also require additional physiological constraints on the phasic SC response parameters as assumption of different phasic SC response parameters will result in too many degrees of freedom in the optimization problem leading to infinitely many solutions.[31], our single channel deconvolution [32], and our concurrent deconvolution, respectively.The subpanels respectively depict recorded SC signal from hand, SC signal from foot, inferred autonomic nervous activity u, stress state w j |J and probability of p j |J with their confidence intervals.The color-coded backgrounds green, light violet, light red, and yellow correspond to rest period, city driving, toll road, and highway, respectively.
Deconvolution of SC signal is a challenging problem as multiple sets of physiological parameters and neural stimuli exist which closely approximate an observed signal.Besides, the smallest level of noise can perturb the solution to a physiologically infeasible point because of the sensitive nature of the system impulse response with respect to the physiological parameters.The proposed formulation is a nonconvex problem and we solve it using a coordinate descent deconvolution approach until convergence to a local minimum.In order account for non-convexity, we use multiple initializations and choose the solution that minimizes the cost function compared to all other solutions that the algorithm finds.To ensure identifiable solution, we use appropriate physiologically plausible constraints similar to [1], [37] on the unknowns.Figure 5 shows the results from previous two algorithms: LedaLab [17] and cvxEDA [19].These two existing methods can solve the inverse problem of finding u using single channel data assuming known physiological parameters.In contrast, the proposed approach can solve for both the physiological system parameters and the inverse problem of finding neural stimuli using multichannel recording.
From our analysis of the auditory stimulation data, we observe a very large phasic response right after a stimulation has been given to the participant.For example, female participant 1, male participant 2, and male participant 3 have significantly distinct phasic SC responses right after the auditory stimuli.Although the phasic SC data from female participant 2 shows multiple responses after one auditory stimulation, the very first response right after an auditory stimulation is usually very prominent.Other comparatively smaller responses indicate that the participant requires more perspiration in order to reduce the body temperature increases due to increased metabolism.Our algorithm successfully detected these small responses as well.In general, the distance between two consecutive phasic responses is more than a few seconds.Therefore, we chose a minimum separation of 1 second between two adjacent peaks in the deconvolution algorithm to obtain a physiologically feasible solution.
SC data can be noisy and small noisy peaks are somewhat comparable to the small insignificant phasic SC responses.To avoid detecting noise peaks as SC phasic response, we used an internal threshold in each iteration GCV-FOCUSS+.Any detected estimated non-zero element of u that is smaller than the threshold is set to zero.In this study, we used 5 as the threshold.This threshold works well for FIGURE 10. White Gaussian structure in the model residual errors of phasic SC DATA of six female participants and six male participants for the recordings corresponding to hand.Each of the panels displays the quantile-quantile plot of the Phasic SC model residual errors; the graph shows that the residual errors are Gaussian.almost all of the participants.However, there are cases where participants have very small phasic SC responses that are very comparable to noise and our algorithm might discard them.For instance, male participant 3 has very small phasic SC responses during the auditory stimulation for both hand and foot.In this case, the phasic SC responses that are comparable to the noise are discarded.Figure 13 in Appendix shows that the simulation results from simulation corresponding to male participant 3 discards one detected pulse that is very close to 5.This is because the addition of noise to simulated data made the smallest phasic SC response comparable to noise.This also explains the low R 2 values in male participant 3.In our future work, we plan to include participant dependent threshold selection to enable us to detect phasic SC responses for cases where the participants have very small ANS activation.In the case of male participant 4, one small peak present between 260 seconds and 280 seconds of foot data is not visible in hand data.This denotes that the peak is due to noise.Our algorithm discards this peak and prevent overfitting to the noise.This is the reason why the R 2 value in male participant 4 foot data is small.
Figure 6 shows how noise can deteriorate estimation accuracy.Figure 6 also shows that with the addition of the more channels, the estimation accuracy of decay time τ d increases.However, the same is not visible for τ r .This might be because τ r is smaller than 1 second which is less than the sampling frequency.It is hard to capture the information of about τ r with 1 Hz sampling frequency.Nevertheless, with higher sampling frequency, estimation accuracy for τ r should also improve with addition to more channels similar to τ d .The stress results obtained in Figure 7 for all the cases are consistent.The result in Figure 7c shows much smoother estimation compared to the result in 7b and 7a.Although stress estimation using the results from single-channel deconvolution obtains a smoother estimate compared to the peak detection one, it is not as good as the concurrent deconvolution.In the first rest period, there is an unwanted peak in the estimated stress state in the heuristic approach.The stress estimation using the result from single-channel deconvolution also detected the unwanted stress state spike.In contrast, the result in Figure 7c shows much-improved stress tracking with no significant stress state spike in the first rest period.A similar spike in the HAI can also be seen in the first rest period for 7b and 7a which is very small in 7c.Simple peak detection algorithm might not provide the accurate timings and amplitudes of the neural stimuli.Moreover, small noise peaks can be captured with the peak detection algorithm where they might not represent actual neural stimuli.

V. CONCLUSION
In this study, we proposed a physiological state-space model for multi-channel SC recordings from different regions of skin.Morover, we prosposed a concurrent deconvolution algorithm for simultaneously collected multi-channel SC data.Analysis with experimental and simulated data showed that algorithm successfully recovers the neural stimuli due to the known auditory stimulation times.Our proposed method and algorithm results in integrating multiple simultaneously collected SC data to recover the ANS stimuli robustly in the presence of noise and different artifacts.Moreover, we applied our approach to concurrently deconvolve simultaneously recorded signals from multiple skin region in real-world driving stress condition.Using the concurrently deconvolved driver's stress data, we were able to achieve a better estimate of stress states than the previous study.The state-space model formulation and deconvolution algorithm successfully recover the stimuli.
As one of our future direction, we plan to use an inference framework to obtain confidence intervals and reduce the time complexity by avoiding large matrix inversions.Furthermore, we plan to include tonic component modelling for multi-channel approach to concurrently separate tonic and phasic components, similar to previous single channel approaches [19], [41].To effectively and robustly estimate psychological stress level, we plan to utilize the inferred ANS activity from multi-channel SC recordings with other physiological signals similar to [42]- [49].Finally, we intend to device appropriate control strategy to reduce the stress level estimated using the inferred ANS activity from a multichannel approach similar to [50]- [53].The panels show the deconvolution result on the simulated data for six male participants, respectively.In each panel, i) the top sub-panel shows the simulated (blue stars) and the estimated (red curve) phasic component corresponding to the middle phalanx of hand; ii) middle panels shows the simulated (blue stars) and estimated (red curve) phasic component corresponding to the the medial planar surface of foot; and iii) bottom sub-panel shows the timings of the simulated ANS activation timings and amplitudes (gray line) and the estimated ANS activation timings and amplitudes (red dashed line).

APPENDIX ADDITIONAL RESULTS
Figure 8 and 9 show deconvolution results from six female and six male participants, respectively.Figure 10 and 11 shows the quantile-quantile plots of the phasic SC model residual errors for the 6 female participants and 6 male participant from both channels.The quantile-quantile plots suggests that the model captures the Phasic SC dynamics, and the phasic SC residual errors have a white Gaussian structure.Figure 12 and 13 show deconvolution results from simulated phasic SC data for six female and six male participants, respectively.

FIGURE 1 .
FIGURE 1. Model Block Diagram.A single neural stimuli signal u(t ) generated by the ANS is responsible for the phasic response in different regions of the skin throughout the body.The block diagram shows the same neural stimuli u(t ) stimulating χ different regions with delay parameters β n , ∀n ∈ {1, 2, 3, • • • , χ}.The attenuation term α n reflects the ratio of the number of sweat glands in the nth region to that of the reference region.
and Q are unitary matrices and the κ i 's are the singular values of D θ P 1 2

FIGURE 2 .
FIGURE 2. Estimated Deconvolution of the Experimental Phasic SC Signals Two Female and Two Male Particiapants: In each of the panels, i) the top sub-panel shows the experimental (red stars) and the estimated (green curve) phasic component corresponding to the middle phalanx of hand; ii) middle sub-panel shows the experimental (blue stars) and estimated (green curve) phasic component corresponding to the the medial planar surface of foot; and iii) bottom sub-panel shows the timings of the auditory stimuli (gray vertical lines) and the estimated ANS activation timings and amplitudes (green vertical lines).

FIGURE 3 .
FIGURE 3. Three Channel Deconvolution on Experimental Data: In each of the panels, i) the first sub-panel shows the experimental (red stars) and the estimated (green curve) phasic component corresponding to the middle phalanx of hand; ii) the next sub-panel shows the experimental (black stars) and estimated (green curve) phasic component corresponding to the thenar/hypothenar of hand; iii) the next sub-panel shows the experimental (blue stars) and estimated (green curve) phasic component corresponding to the medial planar surface of foot; and iv) bottom sub-panel shows the timings of the auditory stimuli (gray vertical lines) and the estimated ANS activation timings and amplitudes (green vertical lines).

FIGURE 6 .
FIGURE 6. Noise vs Accuracy Plot for Rise Time and Decay Time.Left and right sub panels show the how percentage error increases with low SNR.Each data point corresponds to the average percentage error of eight simulated trials.The model parameters for data simulated data used are τ r = 0.75 seconds, τ d = 4 seconds, α 2 = 0.5 and α 3 = 0.3.neural stimuli used for simulation is same as in Figure 5.

FIGURE 7 .
FIGURE 7. Stress State Estimation from Drivers Stress Dataset.All three panels from left to right show estimated stress from the drivers stress dataset with heuristic peak detection[31], our single channel deconvolution[32], and our concurrent deconvolution, respectively.The subpanels respectively depict recorded SC signal from hand, SC signal from foot, inferred autonomic nervous activity u, stress state w j |J and probability of p j |J with their confidence intervals.The color-coded backgrounds green, light violet, light red, and yellow correspond to rest period, city driving, toll road, and highway, respectively.

FIGURE 8 .
FIGURE 8.Estimated Deconvolution of the Experimental Phasic SC Signals Six Female Particiapants: In each of the panels, i) the top sub-panel shows the experimental (red stars) and the estimated (green curve) phasic component corresponding to the middle phalanx of hand; ii) middle sub-panel shows the experimental (blue stars) and estimated (green curve) foot phasic component corresponding to the the medial planar surface of foot; and iii) bottom sub-panel shows the timings of the auditory stimuli (gray vertical lines) and the estimated ANS activation timings and amplitudes (green vertical lines).

FIGURE 9 .
FIGURE 9.Estimated Deconvolution of the Experimental Phasic SC Signals Six Male Particiapants: In each of the panels, i) the top sub-panel shows the experimental (red stars) and the estimated (green curve) phasic component corresponding to the middle phalanx of hand; ii) middle sub-panel shows the experimental (blue stars) and estimated (green curve) foot phasic component corresponding to the the medial planar surface of foot; and iii) bottom sub-panel shows the timings of the auditory stimuli (gray vertical lines) and the estimated ANS activation timings and amplitudes (green vertical lines).

FIGURE 11 .
FIGURE 11.White Gaussian structure in the model residual errors of phasic SC data of six female participants and six male participants for the corresponding to foot.Each of the panels displays the quantile-quantile plot of the Phasic SC model residual errors; the graph shows that the residual errors are Gaussian.

FIGURE 12 .
FIGURE 12.Estimated Deconvolution of the Simulated Phasic SC Signal: The panels show the deconvolution result on the simulated data for six female participants, respectively.In each panel, i) the top sub-panel shows the simulated (blue stars) and the estimated (red curve) phasic component corresponding to the middle phalanx of hand; ii) middle panels shows the simulated (blue stars) and estimated (red curve) phasic component corresponding to the the medial planar surface of foot; and iii) bottom sub-panel shows the timings of the simulated ANS activation timings and amplitudes (gray line) and the estimated ANS activation timings and amplitudes (red dashed line).

FIGURE 13 .
FIGURE13.Estimated Deconvolution of the Simulated Phasic SC Signal: The panels show the deconvolution result on the simulated data for six male participants, respectively.In each panel, i) the top sub-panel shows the simulated (blue stars) and the estimated (red curve) phasic component corresponding to the middle phalanx of hand; ii) middle panels shows the simulated (blue stars) and estimated (red curve) phasic component corresponding to the the medial planar surface of foot; and iii) bottom sub-panel shows the timings of the simulated ANS activation timings and amplitudes (gray line) and the estimated ANS activation timings and amplitudes (red dashed line).
MD. RAFIUL AMIN (S'15) received the B.S. degree in electrical and electronic engineering from the Bangladesh University of Engineering and Technology, in 2016.He is currently pursuing the Ph.D. degree in electrical and computer engineering with the University of Houston, Houston, TX, USA.He worked as a Software Engineer at Samsung R&D Institute Bangladesh, from 2016 to 2017.ROSE T. FAGHIH (S'12-M'15) received the B.S. degree (Hons.)(summa cum laude) in electrical engineering from the University of Maryland, College Park, MD, USA, in 2008, and the S.M. and Ph.D. degrees in electrical engineering and computer science with a minor in mathematics from the Massachusetts Institute of Technology (MIT), Cambridge, MA, USA, in 2010 and 2014, respectively.She completed her Postdoctoral training at the Department of Brain and Cognitive Sciences and the Picower Institute for Learning and Memory, MIT, and the Department of Anesthesia, Critical Care and Pain Medicine, Massachusetts General Hospital.At MIT, she was with the Laboratory for Information and Decision Systems and the MIT-Harvard Neuroscience Statistics Research Laboratory.Her research interests include medical cyber-physical systems as well as control, estimation, and system identification of biomedical systems.She is currently an Assistant Professor with the Department of Electrical and Computer Engineering, University of Houston, where she directs the Computational Medicine Laboratory.Dr. Faghih was a recipient of various awards, including the 2016 IEEE-USA New Face of Engineering Award, the National Science Foundation Graduate Research Fellowship, the MIT Graduate Fellowship, and the University of Maryland's Department of Electrical and Computer Engineering Chair's Award.She has also been inducted into various honor societies, including Phi Kappa Phi, Tau Beta Pi, and Eta Kappa Nu.

TABLE 3 .
Deconvolution errors with our single channel and concurrent deconvolution using simulated data.