A mixed filter algorithm for sympathetic arousal tracking from skin conductance and heart rate measurements in Pavlovian fear conditioning

Pathological fear and anxiety disorders can have debilitating impacts on individual patients and society. The neural circuitry underlying fear learning and extinction has been known to play a crucial role in the development and maintenance of anxiety disorders. Pavlovian conditioning, where a subject learns an association between a biologically-relevant stimulus and a neutral cue, has been instrumental in guiding the development of therapies for treating anxiety disorders. To date, a number of physiological signal responses such as skin conductance, heart rate, electroencephalography and cerebral blood flow have been analyzed in Pavlovian fear conditioning experiments. However, physiological markers are often examined separately to gain insight into the neural processes underlying fear acquisition. We propose a method to track a single brain-related sympathetic arousal state from physiological signal features during fear conditioning. We develop a state-space formulation that probabilistically relates features from skin conductance and heart rate to the unobserved sympathetic arousal state. We use an expectation-maximization framework for state estimation and model parameter recovery. State estimation is performed via Bayesian filtering. We evaluate our model on simulated and experimental data acquired in a trace fear conditioning experiment. Results on simulated data show the ability of our proposed method to estimate an unobserved arousal state and recover model parameters. Results on experimental data are consistent with skin conductance measurements and provide good fits to heartbeats modeled as a binary point process. The ability to track arousal from skin conductance and heart rate within a state-space model is an important precursor to the development of wearable monitors that could aid in patient care. Anxiety and trauma-related disorders are often accompanied by a heightened sympathetic tone and the methods described herein could find clinical applications in remote monitoring for therapeutic purposes.


Mean and variance of the posterior density function p(x k |y k )
We follow an approach similar to [1] in deriving the filter updates. Recall that linear models are assumed to relate sympathetic arousal to the phasic-derived and tonic components.
We take the two noise terms v k and w k to be independent of each other. Consequently, the density functions p(r k |x k ) and p(s k |x k ) conditioned on already having observed x k are independent of each other in the following derivation.
p(x k |y k ) = p(x k |y k−1 )p(y k |x k ) p(y k |y k−1 ) = p(x k |y k−1 )P (m k |x k )p(r k |x k )p(s k |x k )P (n k,1:J |x k ) p(y k |y k−1 ) log(λ k,j ∆)n k,j − λ k,j ∆ + const (5) We take the partial derivative of the logarithm term above and set it to 0 to solve for the mean.
Solving for x k in the equation above provides the filter update for x k|k . We have taken, when calculating the partial derivative. Similarly, the second partial derivative is, The filter update for σ 2 k|k is given by [1], 1.2 Derivation of the M-step updates

Complete data log-likelihood
Taking X K = {x 1 , x 2 , . . . , x K }, the complete data likelihood conditioned on the model parameters Θ is given by, The expected log-likelihood is, (11) Following [2], we take

M-step updates for α and ρ
Let Q 1 denote the term in Q that contains α and ρ.
While it is possible to determine the starting state x 0 as a separate parameter, we follow one of the options in [2,3] and set x 0 = x 1 . This permits some bias at the beginning. Therefore, We take the partial derivatives of Q 1 with respect to α and ρ and set them to 0 to obtain the M-step updates.
The solutions to these simultaneous equations provide α and ρ.

M-step updates for
Let Q 2 denote the term in Q containing γ 0 and γ 1 .
Taking the partial derivatives with respect to γ 0 and γ 1 yields, The solutions to these simultaneous equations provide γ 0 and γ 1 . δ 0 and δ 1 may be obtained similarly from the term in Q containing s k .
1.2.4 M-step updates for σ 2 v and σ 2 w Let Q 3 denote the term in Q containing σ 2 v .
We take the partial derivative with respect to σ 2 v and set it to 0.
The update for σ 2 w may be obtained likewise.

M-step update for σ 2 ε
Let Q 4 denote the term in Q containing σ 2 ε .
We take the partial derivative with respect to σ 2 ε and set it to 0.
We perform a Taylor expansion of the logarithm term around x k|K [1].
Taking the expected value on both sides, Therefore, Now, And similarly, Taking the partial derivative w.r.t. β 0 , And similarly for β 1 we arrive at, 1.2.7 Approximation for the expectation term containing λ k,j Let Q 6 denote the expectation term containing λ k,j .
We perform a Taylor expansion of the summed term around x k|K [1].
Taking the expected value on both sides, Therefore, 2 Experimental data -model parameter estimates The experimental model parameters estimated for each participant are shown in Table 1. Recall that we estimate x k at the E-step and calculate the model parameters at the M-step. Recall also that due to computational complexity we split the estimation of the model parameters related to heart rate (i.e., the θ i 's and the η coefficient) into two parts and calculate them separately. We calculate the θ i 's offline based on maximum likelihood estimation (MLE) and select η based on which value maximized a log-likelihood term. As pointed out in the "Discussion" section of the main text, this separated-out calculation is a limitation of our model (e.g. it can result in numerical issues). The separate estimation of the heart rate parameters may be the reason why η values are small in the final estimates and why larger η values cause convergence issues in the Newton-Raphson solution for the state update x k|k since the MLE estimation of the θ i 's may account for much of the heart rate variability. The value of β 1 also turned out to be negative for two participants (the M-step updates turned out to be negative even after the first iteration). Our algorithm provides two options for calculating β 0 and β 1 and it is possible to select the alternate option which sets β 0 = 1 and calculates β 1 empirically as well if this is to be avoided beforehand. As noted in the main text, a model with a less complex form of the conditional intensity function may enable all the parameters to be recovered at once at the M-step. Lower computational load is also likely to have the benefit of easier deployment onto a wearable platform.