article: Learning to predict threat from environmental cues is a fundamental skill in changing environments. This aversive learning process is exemplified by Pavlovian threat conditioning. Despite a plethora of studies on the neural mechanisms supporting the formation of associations between neutral and aversive events, our computational understanding of this process is fragmented. Importantly, different computational models give rise to different and partly opposing predictions for the trial-by-trial dynamics of learning, for example expressed in the activity of the autonomic nervous system (ANS). Here, we investigate human ANS responses to conditioned stimuli during Pavlovian fear conditioning. To obtain precise, trial-by-trial, single-subject estimates of ANS responses, we build on a statistical framework for psychophysiological modelling. We then consider previously proposed non-probabilistic models, a simple probabilistic model, and non-learning models, as well as different observation functions to link learning models with ANS activity. Across three experiments, and both for skin conductance (SCR) and pupil size responses (PSR), a probabilistic learning model best explains ANS responses. Notably, SCR and PSR reflect different quantities of the same model: SCR track a mixture of expected outcome and uncertainty, while PSR track expected outcome alone. In summary, by combining psychophysiological modelling with computational learning theory, we provide systematic evidence that the formation and maintenance of Pavlovian threat predictions in humans may rely on probabilistic inference and includes estimation of uncertainty. This could inform theories of neural implementation of aversive learning.
  Learning to predict threat from environmental cues is a skill found in many species across the animal kingdom. A laboratory example is Pavlovian threat conditioning (also termed fear conditioning [1]) in which the contingent presentation of predictive cues (conditioned stimuli, CS) and aversive events (unconditioned stimuli, US) engages a process of associative learning [2]. In mammals, including humans, establishing CS/US associations requires synaptic plasticity in basolateral and central amygdala [2–5] and thus relies on a neural circuit distinct from that involved in learning reward associations [6]. Despite progress in elucidating systems-level mechanisms that pre-process and relay CS and US information to the amygdala [7,8], a computational understanding of aversive learning remains incomplete, and available neural data do not fully fit standard learning models [4,9, 10]. Associative learning theory offers a range of computational models, which make specific behavioural predictions. However, many of these models have historically been developed to capture behavioural and neural phenomena in reward learning [11], and it is not known to what extent threat learning follows the same algorithms. Here, we sought to arbitrate between different associative learning models in humans by comparing their trial-by-trial predictions to the measured trajectory of ANS responses, and examine which learning quantities are reflected on different ANS output. We considered models that learn transition probabilities (i. e. parameters) in a known environmental structure, rather than learning the structure itself [10] or its underlying latent causes [12]. Previous models proposed to explain Pavlovian threat learning are derived from classical reinforcement learning (RL) theory and build on Rescorla-Wagner (RW) [13] and Pearce-Hall (PH) [14] rules. Other studies suggested that a combination of these models best captures aversive learning [15,16] (hybrid model, HM). All of these models provide point estimates of future outcomes, by modifying current predictions with a prediction error, weighted by a learning rate that is either fixed (RW) or determined by a fixed proportion of previous prediction errors (PH and HM). While such models can, in theory, be implemented in simple neural architectures [17], at least in the reward domain they have been challenged because they cannot explain some experimental observations such as latent inhibition [18]. However, latent inhibition is easily accommodated by models in which the learning rate is variable and determined by an estimate of the prediction uncertainty [18]. This uncertainty can be implemented in explicit probabilistic computations, for example in hierarchical models [19], or by Kalman filters in the RL framework, which could be classified as implicitly probabilistic and where uncertainty is encoded in a summary statistic, the Kalman gain [18]. Here, we sought to compare previously proposed algorithms of threat learning with a probabilistic account. Notably, probabilistic RL algorithms [18] and general-purpose, hierarchical probabilistic models [19] contain parameters that are difficult to constrain from biological or psychological principles and thus have to be estimated from measured data. We reasoned that this may be difficult given the relatively low signal-to-noise ratio in ANS estimates. However, if threat learning is adaptive and approximates statistical optimality [20], these parameters could be usefully constrained. This is why we relied on a parameter-free probabilistic model constructed from statistical principles, rather than from neural or cognitive considerations. Specifically, from the perspective of an agent that is fully informed about the task structure, and assumes stationary transition probabilities and trial independence, US occurrence follows a Bernoulli process with a US probability parameter. The model then computes the likelihood over this parameter based on the evidence that has been observed so far. The binomial likelihood takes the functional form of a beta distribution. Hence, we implemented a sequentially updated beta-binomial model [21]. The most common way to behaviourally assess human threat learning is by measuring ANS responses, such as skin conductance responses (SCR, mostly sympathetic) [22] or pupil size responses (PSR, sympathetic and parasympathetic) [23]. These are mediated by neural circuits including the central amygdala and thus directly linked to circuits involved in formation of CS/US associations [24,25]. Here, we optimized signal-to-noise ratio with a statistical framework for psychophysiological modelling (PsPM) [26] that exploits the entire time-series of autonomic measurements rather than a low number of arbitrarily selected data features (e. g. peak-to-trough measurements). We have previously shown that this method allows estimating the ANS input causing SCR [27] and PSR [28] with higher precision than standard approaches. Crucially, the observation function that maps an associative learning model onto measured behaviour is not known. Unlike choice data, it has been suggested that ANS outputs may not reflect outcome expectations [15,16]. Furthermore, given their distinct innervation it is possible that different ANS outputs may reflect different aspects of the same learning mechanism. Hence, we additionally sought to clarify this observation function for both measures. We hypothesised that PSR and SCR follow the same underlying learning algorithm but might link to this algorithm by different observation functions. Specifically, SCR are known to habituate during learning, and intraneural stimulations suggested that this habituation does not occur in the effector organ [29], such that either the learning model or the observation function needs to account for this phenomenon.
  To explain trial-by-trial changes in ANS responses, we considered a variety of learning models and different mappings from model onto ANS response (Table 3, Fig 2, S1 Fig). We included a Rescorla-Wagner (RW) model [13] which in our paradigm makes similar predictions as a more general temporal difference rule [32]. A hybrid RW-Pearce-Hall model (HM) [33] which has been used to specifically explain SCR in threat learning, and augments the RW model with a notion of associability [14]. We then considered different observation functions. For the RW model, the only straightforward observation function contains the US expectation. For the hybrid model, the independent variable in the observation function was either associability as previously proposed (HM1) or US expectation (HM2), and for the probabilistic model it was either US expectation (BM) or a combination of US expectation and estimated uncertainty of the environment (BC). As null models, we considered a situation in which no learning occurs over time (i. e. the CS+/CS- difference is time-invariant). We included a probabilistic model that tracks the uncertainty of the environment and does not distinguish CS+/CS- (UN), and a model that distinguishes CS+/CS- but with no change over time (NL). Model evidence was approximated by Bayesian Information Criterion (BIC), which provides a trade-off between model fit and model complexity [34,35]. We first determined the most likely type of observation function for either autonomic measure, irrespective of the underlying model. To this end, we split models into families, according to whether they reflect quantities that are analogous to the expected outcome (3 models: RW, HM2, BM), a notion of previous surprise (also termed associability [15] model HM1), or a combination of expected outcome and its uncertainty (model BC). A family-based model comparison confirmed the intuition that SCR and PSR are mapped onto different observation functions (Fig 1): SCR was best explained by a combination of quantities, indexed by protected exceedance probabilities (p. x. p.) of > 0. 97 for all three experiments (Fig 3A). On the other hand, PSR was best explained by the expected outcome, (p. x. p. > 0. 99, Fig 3A). Within the three models included in the outcome family (RW, HM2 and BM), PSR was best explained by model BM (p. x. p. = 1. 00). We subsequently performed a model comparison across the entire model space (7 models in total, Table 3). Within this space, we examined which combination of learning and observation models explained the trial-by-trial trajectory of estimated SCR and PSR. As in the previous comparison, SCR data were best explained by the combination model BC (Fig 3B p. x. p. of ≥ 0. 99,0. 92 and 0. 96 for Experiments 1–3, respectively). This model implements an ideal observer and maps a mixture of US expectation (prior mean) and uncertainty onto the observation function (Fig 2 for an illustration of model predictions). The latter term—the estimated uncertainty—decreases over time as more instances of the same CS/US contingency are observed (Fig 2 and Table 3). On the other hand, PSR were best explained by model BM, as indexed by p. x. p of 0. 69 (Fig 3B). This model implements the same Bayesian observer as BC, but a different observation function, which corresponds to the US expectation alone (Fig 2 and Table 3). Thus, it appears that pupil size unambiguously tracks expected outcome of a CS, unlike SCR. A crucial difference between previously suggested model HM1 and our model BC is the inclusion of a decay term. This raises a concern that a non-probabilistic model together with a decay term, although not motivated a priori, may explain the data equally well or better. To rule out such a possibility, we asked whether a combination of UN + HM1 would explain the data as well as UN + BM (i. e. BC). We subtracted the estimates of model UN from the SCR data and fitted models BM and HM1 to the residuals. In all three datasets, model BM explained the residuals significantly better than model HM1 (p. x. p. > 0. 99), confirming the result that a probabilistic model better accounts for the data than previously suggested ones. Next, we sought to relate our results to previous studies that only included non probabilistic models into their analysis [15,16]. Among the RL models previously considered, we found highest evidence for HM1, in which SCR reflect the associability term of the Hybrid RW-PH model (Fig 3C, p. x. p. ≥ 0. 99, for experiments 1–3). This result is in keeping with previous work [15,16]. For PSR, the highest evidence was obtained for model RW, which reflects the estimated CS values (p. x. p. of 0. 47, Fig 3C). While it appears that in our and previous studies HM1 best explains SCR among non probabilistic models, a probabilistic model explained our data decisively better, as highlighted above. Notably, our probabilistic model is normative and thus contains free parameters only in the observation function, in contrast to all other models. There is a possibility that the least complex model wins in the quantitative model comparison but does not fit the data very well. To rule out this concern, we computed the variance explained by each model, which provides a quantification of the model fit irrespective of its complexity. For SCR, the highest explained variance across participants was obtained for models HM1 and BC (Fig 4): BC won in experiments 1 and 3 (21% and 18% explained variance) and HM1 won in experiment 2. For PSR, BM and RW explained most variance (both 17%). These results suggest that the probabilistic model fits the data equally well as previously proposed models but with fewer parameters.
  While humans may show overt behaviour during Pavlovian fear conditioning, such as freezing [50], formal attempts to study such phenomena in humans remain limited [51]. This is why we relied on autonomic readouts, the most common way of quantifying fear in humans [2]. In order to analyse ANS activity, we employed a statistical framework for modelling outputs of the peripheral nervous system, such as the SCR and PSR [28,52]. In this framework, estimates of ANS impulse amplitude are obtained using independently validated canonical functions that map ANS firing bursts to peripheral output, with a priori defined constraints on ANS burst timing. Presumably, these estimates more directly reflect central neural states than the noisy peripheral responses themselves. Indeed, ANS amplitude estimates derived from this approach have been shown to better distinguish CS+ and CS- than other approaches, not only for SCR [30,52] and PSR [31] but also for fear-conditioned bradycardia [53], respiratory amplitude responses [54] and fear-potentiated startle [55]. Whether these latter readouts reflect similar quantities as SCR and PSR, and are controlled by the same learning algorithm, remains a topic for future investigation. Our results extend previous fear conditioning studies which were restricted to Rescorla-Wagner and hybrid Rescorla-Wagner-Pearce-Hall models [15,16,56]. Importantly, these studies demonstrated that changes in trial-by-trial SCR can be better explained by an associability term of a RW-PH hybrid model than the US prediction itself [15,16]. When considering non probabilistic RL models only, we could replicate these results, showing that our approach is indeed sensitive and able to capture the same learning mechanisms as previously reported. By using more precise trial-by-trial ANS response estimates and including a probabilistic model, we highlight a different underlying learning mechanism, but conceptually confirm that SCR do not reflect US predictions. In our winning model, the observation function is a combination of the expected US (prior mean) and estimated uncertainty, which depends on the total number of observations during the experiment. Interestingly, there is qualitative evidence that PSR may be informative of past prediction errors [57], something that we did not confirm here. More elaborate experimental designs with specific trial sequences may be required to shed light on this question. Our probabilistic model is normative under stationary outcome probabilities and allows incorporating prior information in the form of beta distributions. This formalism is mathematically elegant but even under an in-built assumption of stationarity it may appear too specific as a candidate for a neural model. Furthermore, there is weak evidence towards a decaying CS+/CS- difference over time, something that requires future experiments to be determined and that none of the models considered here could predict. While our probabilistic model makes the strong assumption of a logarithmic decay of SCR over time, this prediction was not quantitatively corroborated in the current data set. We note that Bayesian models with more generic distributional assumptions such as a Kalman temporal difference filter [18] or an hierarchical Gaussian filter [19] make very similar predictions under the experimental circumstances employed here, such that more sophisticated procedures will be required to disentangle different Bayesian models. Future work may thus elucidate plausible neural implementations that predict the progression of ANS responses highlighted here, which may be distinguished by hemodynamic or electrophysiological activity to formalise the neural implementation of aversive learning.
  The study was performed in accordance with the Declaration of Helsinki, and all procedures for all experiments were approved by the governmental ethics committee (Kantonale Ethikkomission Zurich). All participants gave written informed consent using a form approved by the ethics committee. We conducted three fear conditioning experiments using delay and trace conditioning, and visual as well as auditory CS. A fourth control experiment was used for optimising SCR analysis. Four independent samples of participants (overall N = 102) were recruited from the general and student population (see Table 1 for a detailed sample description). The fourth experiment was identical to experiment 1, and was carried out in a separate group of participants, in order to optimise the estimation of anticipatory sympathetic arousal in an unbiased way. All experiments took place in a dark, soundproof chamber with background illumination provided by the camera and the monitor lights. Participants’ heads were positioned on a chin rest with 70 cm distance from the monitor (Dell P2012H, 20” set to an aspect ratio of 5: 4,60 Hz refresh rate). Skin conductance was recorded from the thenar/hypothenar of participants' left hand, using 8 mm Ag/AgCl cup electrodes (EL258, Biopac Systems Inc. , Goleta, CA, US) and 0. 5% NaCl gel (GEL101, Biopac) [58]. Skin conductance signal was amplified with an SCR coupler/amplifier (V71-23, Coulbourn Instruments, Coulbourn Instruments, Whitehall, PA, US). Data were digitised at 1000 Hz (DI-149, Dataq Instruments, Akron, OH, US), and recorded with Windaq (Dataq Instruments) software. In experiment 3, an EyeLink 1000 System (SR Research, Ottawa, Ontario, Canada) was used for recording pupil diameter and gaze direction with a sampling rate of 500 Hz. Before the experiment we calibrated gaze direction using a 9-point calibration procedure, implemented in the EyeLink 1000 software. All physiological data were analysed using PsPM, version 3. 0. 2 (http: //pspm. sourceforge. net). SCR data were filtered offline with band-pass bidirectional Butterworth filter (cut-off frequencies 0. 0159–5 Hz) and then down-sampled to 10 Hz. We visually inspected the averaged SCR to the US, which normally starts 1–2 s after US onset [59]. Since fear conditioning requires an aversive US, participants who did not show such a US response (defined as a positive peak on the averaged SCR responses exceeding 0. 05 μS relative to a baseline of 2 s before CS onset) were conservatively excluded. This removed 7 participants from experiment 1 and 1 participant from experiment 2. Skin conductance was analysed with a non-linear model (dynamic causal model, DCM) of the anticipatory SCR (aSCR). Note that this approach has no conceptual similarity to the connectivity DCMs used in neuroimaging but uses the same statistical machinery for inversion. This procedure infers activity of the sympathetic nervous system, given changes in the recorded SCR signal [59] and provides trial-by-trial estimates of anticipatory sympathetic responses. For this analysis, we used the default DCM pre-processing and inversion, as implemented in PsPM 3. 0. 2 [30,59]. Timing parameters of the neural input model were initially optimised by using data from an independent experiment, as described below. Prior to all SCR analyses, we optimised the psychophysiological model (PsPM) which specifies the timing of sudomotor inputs into the skin/sweat-gland system with respect to CS. To this end, we used data from a fourth, control experiment, which was not further used for the reinforcement modelling part of the study; hence there is no risk of circularity. Our PsPM describes how sympathetic arousal generates sudomotor impulses, which in turn generate SCR. This model is inverted to yield estimates of sudomotor impulse amplitude. This reduces the influence of non-specific SCR (e. g. spontaneous fluctuations) which are unrelated to the experiment and may overlap with experiment-induced SCR. Because the timing of sudomotor impulses is not known a priori, there are different plausible ways to specify a PsPM; each offers a slightly different balance between model realism and influence of SCR unrelated to the experiment, and thus, a possibly different signal-to-noise ratio. To empirically determine the best PsPM specification for the current experiment, we compared various specifications in their ability to distinguish between CS+ and CS-, something we have previously termed predictive validity [26]. In previous publications we defined a window of sudomotor impulses from CS onset to CS offset and loosely constrained the temporal dispersion of these impulses. This is the most unconstrained definition and based only on the fact that anticipatory arousal must occur after the CS identity is known by the observer, and before the outcome occurs [5,59]. Here, we investigated whether constraining the anticipatory time window improved predictive validity. We fixed the onset of the anticipatory window at CS onset and varied its duration from 0 to 3. 5 s in steps of 0. 5 s. For all models, we considered an additional component reflecting US responses, with fixed onset at 3. 5 s (i. e. US delivery or omission). For each of the resulting models we considered the trial type (CS+/CS-) as dependent variable in a linear regression model, while the estimated responses for each participant, together with participant-specific intercepts were considered as predictors. The F-statistic of this model is identical to the squared t-value from a paired t-test. The Residual Sum of Squares (RSS) of this regression was then converted to Bayesian information criterion (BIC) [60]: BIC=n⋅log (1nRSS) +c Where n is the number of observations and c a complexity constant that is the same for all models. This comparison showed that the model with the lowest BIC value consisted of a fixed component at the onset of the CS with an offset of 0 s (BIC = -76. 74). The BIC of the remaining models ranged from -73. 29 to -59. 19. Therefore, for all subsequent analyses of SCR data, we implemented a DCM with a fixed latency response at CS onset and a fixed latency response at US onset or omission, for each trial. Time-series of pupil size data were extracted from the EyeLink 1000 System, after online parsing for saccade and fixation losses with an in-built algorithm. All data points during eye blinks, head movements, or periods where gaze positions along the x or y axis deviated more than ± 5° visual angle from the fixation, were removed and linearly interpolated [28]. For each participant, we analyzed either the left or right pupil, for whichever there were more data points available. We then estimated the anticipatory input into the pupil with single-trial general linear convolution models (GLMs). Each regressor in the model was formed by convolving a stick function at CS onset for one trial with a synthetic pupil response function previously developed [28,31]. The GLM included one regressor per single trial. We statistically contrasted physiological estimates across trials and participants using linear mixed effects models (LME, package nlme) in R (www. r-project. org; version 3. 2). We considered CS (CS+/-) and Trial as fixed factors and Participant as random, to account for inter-participant variance: ANS˜1+CS*Trial, random=˜1|Participant. LME results were Bonferroni-corrected to account for multiple comparisons across the tested datasets. We assumed that on every trial t, ANS responses, yt, are a linear function of the output of an associative learning model plus a noise term, ϵt: yt=β1⋅zt+β0+ϵt, ϵt∼N (0, σϵ2), where zt is a model variable that links a learning model to the ANS output, while β1 and β0 are model- and participant-specific free parameters. The linear mapping between model output and ANS response is motivated by simplicity, and by the previous observations that in fear conditioning, US probability linear maps onto SCR [61], and that SCR is linearly related to the underlying neural response [62]. We used a Bayesian learning model, which assumes a Bayesian agent that represents threat predictions probabilistically [63,64]. Specifically, a model of the discrete US prediction (Bernoulli probability) with parameter θCS for each of the two CS is sequentially updated based on new sensory information and prior beliefs, according to Bayes’ rule: pt (θ|ut) =p (ut|θ) ⋅pt−1 (θ|ut−1) p (ut) where ut is the presence or absence of US input at trial t. In our case, p (ut | θ), the likelihood function for each individual trial, follows a Bernoulli distribution, as trials can have two possible outcomes (US or no US): p (ut|θ) =θut (1−θ) 1−ut, where ut ∈ {0,1}. The prior probability distribution before the first trial, p0 (θ | ut−1) was assumed to follow a Beta distribution, which is conjugate to the Bernouilli, and which has a probability density function: p (θ|ut−1) =θαt−1 (1−θ) βt−1Β (αt, βt), where B (αt, βt) is a Beta function with parameters αt and βt. The resulting posterior distribution is thus also a Beta distribution, whose parameters are updated according to: αt = αt−1 + ut−1 and βt = βt−1 − ut−1 + 1, where ut = 1 if a US occurred, and ut = 0 otherwise. We assumed uninformative initial priors for the Bayesian model with α0 = β0 = 1, such that the model essentially reflects an ideal Bayesian learner. Although other probabilistic models with more realistic distributional assumptions can account for learning in a Bayesian framework [18,65], they will generate very similar predictions when the prior distribution is uninformative. Hence, we used the Beta-Binomial model for the sake of simplicity: the resulting conjugate Beta posterior distribution has a support in (0,1) and can readily account for the probability of a US to occur. This model allows for different observation functions (Table 3). We considered (a) the mean of the prior distribution (BM), and (b) a combination of (a) with the uncertainty of the environment (BC). a. Prior mean (BM) This response function maps E[θ], the mean of the prior distribution, onto ANS responses: E[θ]=∫θθ⋅p (θ|ut−1) dθ. For the case of a beta distribution the model output, mapped on the independent variable zt, corresponds to: zt=E[θ]=αt−1αt−1+βt−1. b. Combined prior mean and prior uncertainty model (BC) The logarithm of total amount of observations drawn by participants has been previously used to quantify the uncertainty in a changing environment [38]. In our case, where the environment is static and contingencies do not change over the time-course of the experiment, this term is related to the prior uncertainty, or in other words, to the sharpness of the prior distribution. This quantity was computed for each CS separately. A wide distribution corresponds to large values of vt and would allow for large changes in the mean of this distribution from one trial to the next [38]. We considered a model whose independent variable, zt, comprises a mixture of the prior mean and uncertainty: zt=vt+E[θ]. We included into our analysis a range of classical, non-probabilistic reinforcement learning models in which point estimates of US predictions are updated according to some learning rule that involves signed or unsigned prediction errors. a. Rescorla-Wagner rule (RW) The RW model [13] specifies that US prediction (' associative strength' ) is updated according to a signed prediction error signal (i. e. the difference between the prediction on a given trial, xt and the observed sensory input, ut), weighted by a fixed learning rate, η which is a subject-specific free parameter across both CS, such that such that 0<η<1: xt=xt−1+η⋅ (ut−1−xt−1). In our case, where there is only one type of US, which is either present or absent, ut can take two values: ut = 1, if a US was presented on a given trial and ut = 0 otherwise. We assumed that participants had no specific expectations about the two stimuli and thus set the starting values x0 = 0. 5 for both CS+/CS-. For the RW model we used only one observation function (Table 3): zt=xt. b. Hybrid RW/PH model (HM) The Rescorla-Wagner and Pearce-Hall models can be combined into a hybrid model (HM), which assumes signed prediction errors (like a RW rule), and an extra associability term, ηt [15,16,66]. This quantity, ηt, reflects a dynamic learning rate, which is updated over the course of the experiment according to: ηt=k⋅|xt−1−ut−1|+ (1−k) ⋅ηt−1, xt=xt−1+ηt−1⋅ (ut−1−xt−1), where k is a participant-specific positive free scaling parameter common to both CS, reflecting how fast each participant’s predictions are updated, with 0<k<1. As for the RW model, we assumed that participants had no specific expectations about the two stimuli and thus set the initial values at x0 = 0. 5 for both CS+/CS-. HM could have two different observation functions to ANS responses: the current associability (HM1), and the current US prediction (HM2) Hence, the two possible observation functions of the hybrid model are (Table 3): zt=ηt, zt=xt. We additionally considered two null models. The first of them, model VO, only reflects the total amount of observations drawn by participants: zt=−ln (αt+βt). This model could explain effects of habituation, in the absence of any learning. Finally, we included a null model, which assumes that no habituation takes place over the course of the experiment, and that learning takes place prior to the experiment (i. e. the CS+/CS- difference is constant over time): zt={1,0}. Models were optimised for each participant separately, by minimising the residual sum of squares (RSS) between each model’s predictions and trial-by-trial estimates of anticipatory neural responses: RSSMM=∑t=1T (β1⋅zt+β0−yt) 2, where T is the total number of non-reinforced CS+ and CS- for each experiment. All trials were used for generating model predictions, but only trials which were not paired with a US were taken into account when computing the RSS, in order to avoid contamination by an evoked response to the US, an approach similar to previous studies [15]. The free parameters of each model and two regression parameters of the observation function were fitted using an interior point search algorithm as implemented in the Matlab function fmincon, using the RSS as objective function. Model evidence was approximated by Bayesian Information Criterion (BIC), which was calculated for each model and participant [34,35]: BICMM=p⋅ln (T) +T⋅lnRSSMMT, where p refers to the total number of parameters for each model and T to the total number of observations (i. e. single trials unpaired with a US). BIC approximates the true model evidence and provides a compromise between model fit and model complexity, indicated by the total number of free parameters in each model, p. For all models, p ≥ 2, since at least two parameters were estimated for the linear mapping between transfer functions, g, and the psychophysiological data. Model selection was performed at the group level, using random effects analysis (RFX), [67,68], which treats models as random effects in the population. RFX was based on protected exceedance probabilities (p. x. p.), which quantify the probability of a model to be more likely than any other, give the group data, using SMP12 (http: //www. fil. ion. ucl. ac. uk/spm/software/spm/). RFX originally accounts for the possibility that different individuals use different models. This is unlikely to be the case in an evolutionarily conserved mechanism like threat learning, but RFX also provides good protection against outliers which is why we used it here. To examine the a priori model selection error of the given experimental design and model space, we simulated data based on an already published dataset with an independent and larger sample of participants in the same experimental design [37]. We fitted all of our models to skin conductance data obtained during the acquisition phase from the placebo participant group (N = 38) of this study. This yielded, for each model, an empirical distribution of participant-wise parameters. Then, we simulated 256 experiments with 20 participants each (i. e. a sample size similar to that of the individual experiments in the manuscript), by bootstrapping from the estimated distribution of parameter values. To the simulated datasets we fitted all the models and computed model recovery rates using BIC as the model selection criterion.
summary: Using environmental cues to predict threat is a crucial skill, encountered in many species. A laboratory model to study associating predictive cues with aversive events is Pavlovian fear conditioning. A computational understanding of this process at a systems-level is lacking. Here, we investigate which computational learning model best predicts activity of the human autonomic nervous system during fear conditioning. We show, across three data sets and two autonomic readouts, that a probabilistic learning model explains the data decisively better than previously proposed models. We suggest that humans learn to predict threat by maintaining and constantly updating a probabilistic model of the environment through Bayes’ rule.
keywords: learning
  applied mathematics
  experimental design
  conditioned response
  social sciences
  mathematical models
  neuroscience
  learning and memory
  simulation and modeling
  algorithms
  research design
  cognitive psychology
  mathematics
  membrane electrophysiology
  bioassays and physiological analysis
  research and analysis methods
  behavior
  mathematical and statistical techniques
  electrophysiological techniques
  fear conditioning
  behavioral conditioning
  psychology
  electrode recording
  biology and life sciences
  physical sciences
  cognitive science
title: Human Pavlovian fear conditioning conforms to probabilistic learning