Whole brain correlates of individual differences in skin conductance responses during discriminative fear conditioning to social cues

Understanding the neural basis for individual differences in the skin conductance response (SCR) during discriminative fear conditioning may inform on our understanding of autonomic regulation in fear-related psychopathology. Previous region-of-interest (ROI) analyses have implicated the amygdala in regulating conditioned SCR, but whole brain analyses are lacking. This study examined correlations between individual differences in SCR during discriminative fear conditioning to social stimuli and neural activity throughout the brain, by using data from a large functional magnetic resonance imaging study of twins (N = 285 individuals). Results show that conditioned SCR correlates with activity in the dorsal anterior cingulate cortex/anterior midcingulate cortex, anterior insula, bilateral temporoparietal junction, right frontal operculum, bilateral dorsal premotor cortex, right superior parietal lobe, and midbrain. A ROI analysis additionally showed a positive correlation between amygdala activity and conditioned SCR in line with previous reports. We suggest that the observed whole brain correlates of SCR belong to a large-scale midcingulo-insular network related to salience detection and autonomic-interoceptive processing. Altered activity within this network may underlie individual differences in conditioned SCR and autonomic aspects of psychopathology.


Introduction
Discriminative fear conditioning is one of the most common methods of studying fear learning in laboratory settings (LeDoux, 2012;Davis andWhalen, 2001 andGrillon, 2002). It refers to the pairing of an initially neutral stimulus (CS+) with an aversive unconditioned stimulus (US), thereby imbuing the CS+ with the capacity to elicit increased autonomic responses relative to an unpaired neutral stimulus (CS-). Acquisition of conditioned fear forms an important part of theories of how anxiety and stress disorders develop from experiencing aversive events (Mineka and Oehlberg, 2008;Bouton et al., 2001;Craske et al., 2014). A more precise understanding of why some individuals are more prone than others to developing anxiety and stress disorders may therefore be achieved by studying individual differences in fear conditioning. Ultimately, such knowledge could inform the development of more effective treatments for subgroups of patients whose responses to fear conditioning are related to the disorder (see e.g. Insel, 2014). Consequently, it has been argued that fear conditioning research should increasingly focus on understanding the sources of individual differences. One way that individuals vary during fear conditioning is in their autonomic responses (Lonsdorf and Merz, 2017).
The most commonly used autonomic measure of conditioned fear is the differential Skin Conductance Response (SCR; , a label referring to the changes in skin conductance induced by sympathetic activation (Dawson et al., 2007;Wallin, 1981). A plethora of studies have investigated SCR in fear conditioning to determine its associations with psychiatric disorders. A metaanalysis of studies in patients with anxiety disorders found that SCR to the control cue (CS-) during the acquisition of conditioned fear was elevated in patients relative to controls (Duits et al., 2015). Duits et al. have proposed that increased SCR to the CS-was an effect of fear generalization from the threat cue to the control cue, or of reduced inhibition of threat responses to the control cue. Also, in individuals with schizophrenia, an increased SCR to the CS-has been reported across four independent fear conditioning studies (Tuominen et al., 2022). In OCD, however, a systematic review found associations to be less clear when considering SCR during the acquisition of conditioned fear and point to stronger OCD-related differences in SCR during the extinction phase of conditioned fear (Cooper and Dunsmoor, 2021). Taken together, these studies demonstrate increased SCR during fear conditioning to the control cue in patients with anxiety disorders and schizophrenia, relative to control participants, which gives credence to the notion that fear conditioning is an experimental model that can inform research on autonomic regulation in these disorders (Lonsdorf and Merz, 2017).
Previous studies of the neural correlates of individual differences in conditioned SCR, generally defined as the difference in average SCR score between CS+ and CS-presentations during acquisition (see Lonsdorf and Merz, 2017, for a discussion of definitions), have focused on either one or a few brain regions, using region of interest analyses. Many of these studies have found positive correlations with neural responses in the amygdala (LaBar et al., 1998;Phelps et al., 2004;Dunsmoor et al., 2011;Petrovic et al., 2008;MacNamara et al., 2015;Marin et al., 2020). Neuroimaging studies of within-subject variation in conditioned SCR have also generally found positive correlations to amygdala responses (Cheng et al., 2003;Knight et al., 2005;Cheng et al., 2006), although exceptions exist (Sjouwerman et al., 2020;Savage et al., 2021). The findings from studies that report a positive relationship between SCR and amygdala activity are in line with the general understanding of fear conditioning from animal models, where a neural circuitry centered on the amygdala is responsible for the acquisition of conditioned fear responses (LeDoux, 2000;Davis, 2000). They also complement those human lesion studies demonstrating either diminished or absent conditioned SCR following amygdala damage (LaBar et al., 1995;Bechara et al., 1995), although not all studies have found such an effect (Ahs et al., 2010; for a review see Ojala and Bach, 2020). Further, the involvement of the amygdala in human fear conditioning has been questioned based on the results of a meta-analysis of fMRI studies investigating fear conditioning (Fullana et al., 2016) and based on studies showing unexpected, increased amygdala responses to the CS-compared to the CS+ (see e.g. Visser et al., 2021). Such results could arise from distributed representations of the CS+ and CS-in the amygdala (Bach et al., 2011;Reijmers et al., 2007) or from a need for larger sample sizes to detect differential responses in the amygdala. Speaking to the latter idea, two independent studies, each including hundreds of participants, have recently reported increased CS+, relative to CS-, activation in the amygdala (Kastrati et al., 2022;Wen et al., 2022). Amygdala activation to the CS+ was primarily detected during the first trials of acquisition, whereas CS-activity was larger in the end of acquisition (Wen et al., 2022). The results of these two large and independent neuroimaging studies, together with the fairly consistent findings of correlated individual differences in conditioned SCR and amygdala activation (LaBar et al., 1998;Phelps et al., 2004;Dunsmoor et al., 2011;Petrovic et al., 2008;MacNamara et al., 2015;Marin et al., 2020), support the hypothesis that amygdala activation should be positively correlated with SCR.
A limiting factor of previous studies of the neural correlates of individual differences in conditioned SCR is that they generally have been based on small sample sizes (N ≤ 27) and have reported results from region of interest analyses (ROIs). Notably, most previous studies have sample sizes that fall below the minimum guidelines for correlation analysis in fMRI research (Yarkoni, 2009;Yarkoni and Braver, 2010), which recommend a sample size of at least N = 40 for determination of inter-individual correlations. Studies that do not comply with this minimum recommendation are highly susceptible to type II errors (i.e. not detecting real effects, even within ROIs) as well as gross inflation of reported effect sizes (Yarkoni, 2009). Furthermore, some of the reported studies also use non-standard statistical procedures, such as reporting results from uncorrected whole brain analyses (Petrovic et al., 2008), or initially use an uncorrected whole brain analysis followed up by stringent FWE-correction within masks only targeting regions implicated in the prior uncorrected analysis (Dunsmoor et al., 2011), which raises the risk for type I error (see e.g. Poldrack et al., 2017) and could reduce the reliability of the results. Only two previous studies that investigated the association between individual differences in conditioned SCR and neural activity have met the minimum requirements for sample size suggested by Yarkoni andBraver, 2010: MacNamara et al., 2015 (N = 49) and Marin et al., 2020 (N = 60). However, both of these studies have still reported results based on uncorrected statistics in pre-defined ROIs (See Table 1 for a comparison of previous studies investigating neural correlates of SCR).
The previous focus on ROIs excludes activations in many parts of the brain that are potentially important for explaining individual differences in conditioned SCR. Fear conditioning is known to activate a large set of cortical, subcortical, and brainstem areas other than the ROIs that have so far been investigated for their correlation with SCR (Fullana et al., 2016). In their meta-analysis, Fullana et al., 2016 proposed that neural regions consistently activated during fear conditioning collectively constitute a large-scale neural network, centered on the dACC and anterior insula, that represents autonomic-interoceptive processing in response to conditioned stimuli (Fullana et al., 2016; based on findings by e.g. Cameron, 2009;Craig, 2009;Critchley and Harrison, 2013;Medford and Critchley, 2010). Based on this proposal, one would expect individual differences in autonomic conditioned responding, such as those measured by SCR (Dawson et al., 2007), to correlate with conditioningrelated activity within this broader network.
Because previous studies only studied correlations to SCR in a handful of brain regions and in relatively small samples, the primary aim of this study was to investigate the whole brain correlations of individual differences in conditioned SCR by analyzing data from a large twin sample performing a fear conditioning task (N = 285 individuals). Similar to previous studies (e.g. Phelps et al., 2004;Dunsmoor et al., 2011;MacNamara et al., 2015), the present study used differential SCR as a between-subjects regressor of CS+ > CS-BOLD activation. Previous studies have not found wholebrain correlations surviving correction for multiple comparisons. In the current study, the sample size of N = 285 individuals provided sufficient statistical power to detect correlations of medium effect size and above (r > .251) throughout the whole brain (see section 4.3.3 in the Materials and methods section for details), thereby substantially improving upon the power of previous whole brain analyses of smaller sample sizes. A family-wise error (FWE) corrected alpha level of α = .05 was used to constrain the risk of type I errors to an acceptable level (see e.g. Poldrack et al., 2017), which was not done in previous studies. Based on findings from previous studies of individual differences in conditioned SCR (LaBar et al., 1998;Phelps et al., 2004;Dunsmoor et al., 2011;Petrovic et al., 2008;MacNamara et al., 2015;Marin et al., 2020), we hypothesized a positive correlation between SCR and amygdala activation. Again, our sample size yielded substantially improved power in comparison to previous studies while securing an alpha level of α = .05 (see section 4.3.3 in the Materials and methods section for details). Finally, a last aim of the present study was to determine whether areas whose activity explained significant individual variation in conditioned SCR did so independently of one other. A comparison of the relative contributions of different brain areas to conditioned SCR has the potential to elucidate separate neural pathways mediating individual differences in conditioned autonomic responding. Ultimately, such findings may aid the understanding of autonomic regulation in pathological fear and anxiety.

Brain responses
The effects of fear conditioning on fMRI responses during habituation and acquisition (CS+ > CS-) are described in Appendix 6. We found no differences in neural responses to the CS+ compared to the CS-during habituation. During acquisition, the pattern of activation to the CS+ relative to the CS-was very similar to the pattern reported in the meta-analysis by Fullana et al., 2016 and included large parts of the striatum, the insula, midline areas of the cingulum, lateral temporal cortex, parietal cortex, and the supplementary motor areas. Of note, the whole brain analysis also revealed greater activation to the CS+ than to the CS-bilaterally in the amygdala.
Correlation between SCR difference scores and brain responses during fear conditioning Whole brain analysis Fear conditioning-related brain responses (CS+ > CS-) were correlated to SCR difference scores in the dorsal anterior cingulate cortex/anterior midcingulate cortex, right anterior insula, right inferior frontal gyrus/frontal operculum, bilateral temporoparietal junction/superior temporal gyrus, right superior parietal lobe/postcentral gyrus, bilateral superior frontal gyri/dorsal premotor cortex, and a right-lateralized midbrain region in areas consistent with periaqueductal gray and reticular formation. For a summary of results, see Figure 1 and Table 2. In order to ensure the reliability of our findings and facilitate comparison with prior studies, we performed the same analysis using average square root transformed raw value SCRs instead of Z transformed SCRs. We also performed the same analysis without any participant exclusion (see the Materials and methods section regarding participant exclusion). Both of these analyses resulted in a very similar pattern of correlations (see Appendix 3 and 4), showing that the results seem robust to different types of SCR normalization and to different choices regarding participant exclusion. In addition, we also repeated our analysis using an extended SCR response window as well as controlling for shock expectancy and genetic influence, again with similar results (see Appendix 7-table 1). Finally, as pointed out by a reviewer, the Ledalab software package (v 3.4.9; Benedek and Kaernbach, 2010) presently used for SCR scoring has been shown to yield no better and sometimes worse results than standard peak scoring (Bach, 2014). For this reason, we also repeated our whole brain analysis using the PsPM software package (v 5.1.1) (Bach and Friston, 2013). This analysis once again implicated the same set of regions. See Appendix 7 for results and Appendix 8 for additional information on how this analysis was conducted (complementary analysis).

Correlation between individual differences in conditioned SCR and amygdala activation
To test the hypothesis of a correlation between individual differences in conditioned SCR and amygdala activation, we performed an ROI analysis. This analysis demonstrated significant correlations bilaterally in the amygdala (right peak MNI coordinates: 20, -2, -14; cluster size = 38 voxels; t = 3.68; left peak MNI coordinates: −22, -2, -16; cluster size = 8 voxels; t = 3.08). See Figure 2. Figure 1. Correlation between individual differences in conditioned SCR and whole brain responses during fear conditioning, obtained using individual SCR scores (Z transformed average CS+ minus CS-SCR) as a second level, between-subjects regressor of the average CS+ > CS-BOLD activation in SPM12 (Wellcome Department of Cognitive Neurology, University College, London) software. The sample consisted of 285 participants who passed the following exclusion criteria: pregnancy, inability to lie still for a 1 hr duration, intolerance of tight confinements, ongoing psychological treatment, metal objects in the body (due to surgery, fragmentation, etc.), current alcohol or drug-related problems, use of psychotropic medications, unsuccessful recording of skin conductance responses, loss of brain imaging data due to excessive head movement, and participant failure to comply with task instruction regarding button press in at least 80% of trials. (A) Activation map of key implicated neural regions. Color-coded t values ranged from t = 3 to t = 6. The statistical image was thresholded at P < 0.05 FWE-corrected and displayed on an anatomical brain template. (B) Scatter plots depicting correlation between SCR difference scores and eigenvariates from significant whole brain clusters in the dorsal anterior cingulate cortex (upper panel) and the temporoparietal junction (lower panel). R = right dACC = dorsal anterior cingulate cortex. TPJ = temporoparietal junction. IFG = inferior frontal gyrus. PAG/RF = periaqueductal gray/reticular formation. AI = anterior insula.
The online version of this article includes the following source data for figure 1: Source data 1. Variable data used to produce Figure 1B, Figure 2B and statistical analyses reported in the section 'Relative contribution of neurofunctional correlates to individual differences in SCR', as well as Appendix 1-figure 1, Appendix 2-figure 1 and Appendix 5-table 1.

Relative contribution of neurofunctional correlates to individual differences in SCR
In order to examine the independent and/or shared contributions of neural responses to SCR, we extracted eigenvariates of contrast values from all significant clusters in the whole brain analysis. Note. MNI coordinates and t values represent significant peak voxels of each cluster. Statistical significance was calculated using t tests implemented within the SPM software with an FWE corrected alpha level of α = .05.

Figure 2.
Correlations between individual differences in conditioned SCR and amygdala activation, obtained using individual SCR scores (Z transformed average CS+ minus CS-SCR) as a second level between-subjects regressor of the average CS+ > CS-BOLD activation in SPM12 (Wellcome Department of Cognitive Neurology, University College, London) software. The sample consisted of 285 participants who passed the following exclusion criteria: pregnancy, inability to lie still for a 1 hr duration, intolerance of tight confinements, ongoing psychological treatment, metal objects in the body (due to surgery, fragmentation, etc.), current alcohol or drug-related problems, use of psychotropic medications, unsuccessful recording of skin conductance responses, loss of brain imaging data due to excessive head movement, and participant failure to comply with task instruction regarding button press in at least 80% of trials. (A) Activation map depicting significant activation on coronal section at MNI Y-coordinate = -2. Color-coded t values range from t = 2.0 to t = 4.0. The statistical image was thresholded at p < 0.05 FWE-corrected. (B) Scatter plot depicting correlation between SCR difference scores and eigenvariates from the significant right amygdala cluster within the amygdala ROI. For source data to (B), see Figure 1-source data 1.
Eigenvariates were then entered as regressors in a hierarchical regression analysis. Together, regional eigenvariates demonstrated a significant correlation with conditioned SCR with a moderate effect size (F(1, 284) = 4.82; r = .37; r 2 = .14; p < 0.001). No region contributed unique, statistically significant variance (see Appendix 5). Source data for all hierarchical regression analyses can be found in Figure  1-source data 1. Individual differences in SCR difference scores could be associated with individual differences in SCR to both the CS+ and the CS-. Therefore, we wanted to test whether SCR to the CS+, and not the CS-, was the reason for the observed correlation between SCR difference scores and eigenvariates. To this end, we correlated the extracted eigenvariates from regions that were correlated to SCR difference scores in the whole brain analysis with average raw value SCRs to the CS+ and CS-, separately. Results demonstrated significant correlations between all regional BOLD eigenvariates and CS+ SCR (p-values ≤ 0.001), except for right superior parietal lobe (p = 0.002), left superior frontal gyrus (p = 0.002), right amygdala (p = 0.006), and left amygdala (p = 0.029), where correlations did not survive the Bonferroni-corrected threshold (α = 0.00151; see the Materials and methods section for details regarding the Bonferroni-correction). However, no correlations between SCR to the CS-and extracted eigenvariates were significant (p-values > 0.00151), and all extracted values were significantly more correlated to CS+ SCR than CS-SCR (p-values < 0.00151) except for the right (p = 0.018) and left amygdala (p = 0.011). This indicated that neural correlations to differential SCR were mainly explained by increased responding to the CS+. Source data for all correlation analyses involving the CS+ and CS-separately can be found in Figure 1-source data 1.

Discussion
Individual differences in SCR during discriminative fear conditioning are common (Lonsdorf and Merz, 2017) and have been associated with psychopathology (Duits et al., 2015;Nees et al., 2015;Lonsdorf and Merz, 2017). In this study, we examined the whole brain correlates of conditioned SCR in a large sample of twins (N = 285) during discriminative fear conditioning with concomitant SCR and fMRI recordings. As expected, we found a correlation between individual differences in conditioned SCR and amygdala activity in line with previous reports (LaBar et al., 1998;Phelps et al., 2004;Petrovic et al., 2008;Dunsmoor et al., 2011;MacNamara et al., 2015;Marin et al., 2020). Our analysis also implicated the dACC and insula, two regions which have previously been found to show increased responses in high vs. low conditioners (Marin et al., 2020). Importantly, we also identified correlations in novel regions including the bilateral temporoparietal junction/superior temporal gyri, right frontal operculum, bilateral dorsal premotor cortex, right superior parietal lobe, and midbrain. All correlations between brain responses and conditioned SCR were positive, meaning that individuals who respond more strongly to the CS+ relative to the CS-on the physiological level (SCR) also showed greater neural activation to the CS+ relative to the CS-. Furthermore, all regional activations demonstrated a stronger correlation with SCR to the CS+ alone compared to the CS-alone, although a few correlations did not survive Bonferroni correction. This indicated that neural activity was primarily associated with heightened, conditioned SCR to the CS+ (i.e. as opposed to inhibited SCR to the CS-).
An important research question is whether the neural network associated with individual differences in conditioned SCR is embedded in the network of regions that is generally activated during fear conditioning. The whole brain correlates of SCR found in the present study belong to the set of regions consistently activated during human neuroimaging studies of fear conditioning (Fullana et al., 2016). This indicates that these regions not only respond to the CS+ relative to the CS-in general, but that the magnitude of this activation also co-varies with individual differences in the magnitude of conditioned responding indexed by SCR. This is consistent with the proposal by Fullana et al., 2016; based on findings by e.g. Cameron, 2009;Craig, 2009;Critchley and Harrison, 2013;Medford and Critchley, 2010 that these regions, especially the dACC and anterior insula, are part of a large-scale neural network regulating autonomic responding (SCR). Reasonably, increased autonomic activation would correlate with larger responses in neural regions regulating autonomic processing. To better understand the potentially unique contributions of different brain regions to SCR, we performed a hierarchical regression analysis. Results suggested that cortical areas together with midbrain regions contributed to individual differences in conditioned SCR. No unique contribution from any of the regions could be proven. These results are consistent with the idea that the regions identified in the whole brain analysis form part of one functional network, rather than separate independent networks. Indeed, the functional connectivity and network-structure between the regions reported in the present study have been examined by several research groups previously (see e.g. Uddin et al., 2019, for a review). These research groups have proposed a more general function for this network, beyond autonomic-interoceptive processing and autonomic regulation, in detecting and preparing responses to salient events across homeostatic, affective, and cognitive domains (see e.g. Seeley et al., 2007;Menon and Uddin, 2010;Uddin, 2015;Menon, 2015;Uddin et al., 2019). Uddin et al., 2019 refer to this network as the 'midcingulo-insular network' to reflect the anatomy of the network, rather than using functional labels that can be dependent on context (e.g. 'salience network', Menon, 2015;'ventral attention network', Corbetta et al., 2008). Within this network, the dACC/aMCC (see Vogt, 2016, for a discussion regarding the naming of this region) and anterior insula, constitute the major input-output nodes (Uddin et al., 2019;Menon, 2015). While, the insula is thought to integrate cognitive, affective, interoceptive, and homeostatic information, the dACC is believed to represent this summarized information in order to determine autonomic, behavioral, and cognitive responding (Menon, 2015;Medford and Critchley, 2010). Efferent autonomic output from the dACC is proposed to be mediated by the periaqueductal gray (Menon, 2015), which was another region identified in our whole brain analysis (labeling consistent with a previous definition of the periaqueductal gray, see Linnman et al., 2012). Notably, the midcingulo-insular network also includes the bilateral temporoparietal junction/superior temporal gyri (Uddin et al., 2019; see e.g. Yeo et al., 2011;Seeley et al., 2007;Bzdok et al., 2013), whose activity we found to be correlated with SCR. The right frontal operculum, which was another region we found to be positively correlated to SCR, is part of what has previously been called the 'cingulo-opercular network' or the 'ventral attention network' (Uddin et al., 2019; see e.g. Corbetta et al., 2008). More specifically, the right lateralized temporoparietal junction and frontal operculum appear to be recruited particularly during exogenous salience detection (Uddin et al., 2019), which fits well with the stimulus-driven salience detection in the context of fear conditioning. Finally, the amygdala has also been proposed to constitute a major subcortical node within this network (Menon, 2015;Uddin et al., 2019). Thus, most of the regions whose activity we found to correlate with SCR belong to this functional network, consistent with its role in regulating autonomic responses, the only exception being the bilateral dorsal premotor cortex. However, the premotor cortex is known to be a major cortical elicitor of the SCR (Dawson et al., 2007;Boucsein, 2012), making it a plausible link between midcingulo-insular activity and peripheral skin conductance. Consequently, based on our findings and on previous evidence-based theory, we propose that individual differences in conditioned SCR may originate from activity within the midcingulo-insular network (Uddin et al., 2019) in conjunction with the dorsal premotor cortex.
It has been suggested that understanding sources of individual differences in fear conditioning may uncover individual risk and resilience factors with respect to fear and anxiety that may ultimately aid the understanding and treatment of fear-related psychopathology (Lonsdorf and Merz, 2017). Our proposal that individual differences in conditioned SCR co-vary with activity in a large-scale neural network substantiating autonomic-interoceptive processing and salience detection may highlight such a source. Indeed, overactivity within the midcingulo-insular network has previously been suggested to underlie pathological anxiety (Menon, 2011;Menon, 2015), as patients with anxiety disorders show altered functional connectivity within this network (Peterson et al., 2014) as well as hyperactivity within the anterior insula (Paulus and Stein, 2006;Stein et al., 2007) and amygdala (Etkin and Wager, 2007), which are important nodes of the network. Based on these findings, it has been suggested that anxiety disorders and neuroticism may result from excessive processing of emotion-related salient cues (Menon, 2011) and/or alterations in autonomic interoceptive-autonomic processing (e.g. Paulus and Stein, 2006;Medford and Critchley, 2010). As anxiety disorders have also been associated with altered discriminative fear learning (Nees et al., 2015), and as our results indicate that individual differences in discriminative fear learning covary with midcingulo-insular activity, our results are consistent with this anxiety model. Specifically, in reviewing the relationship between discriminative fear learning and anxiety disorders, Nees et al., 2015 note that while discriminative fear learning does not appear to be impaired across all anxiety disorders and across all stimuli types (see also Duits et al., 2015, for a meta-analysis confirming this observation), studies comparing anxiety disorder patients to controls with regard to disorder-specific stimuli have found increased discriminatory fear learning in patients (in specific phobia, see Schweckendiek et al., 2011; in social phobia, see Lissek et al., 2008; in PTSD, see Wessa and Flor, 2007), suggesting that discriminatory fear learning may be a mechanism underlying these disorders. The present finding, that individual differences in discriminatory fear learning covary with differences in midcingulo-insular activity, suggests that alterations in such neural activity may also be contributing to, or resulting from, anxiety, consistent with the anxiety model proposed by Menon, 2011. However, it should be noted that the present study only used SCR as an outcome measure of discriminative fear learning while the previously cited studies considered other outcome measures as well (e.g. fear potentiated startle in Lissek et al., 2008; subjective ratings of valence and arousal in Wessa and Flor, 2007), thus somewhat limiting the generalizability of our findings (for further details regarding methodological differences and similarities we refer to the review by Nees et al., 2015). We recommend that future studies continue exploring the midcingulo-insular network and its relationship to fear and anxiety disorders.
One limitation of the present study is the use of social stimuli as CSs. While the general CS+> CS-BOLD contrast analyzed in the present study largely demonstrated a pattern of activation typical to fear conditioning (Fullana et al., 2016; see Appendix 6), the potential influence of social threat processing cannot be entirely ruled out. Therefore, we recommend reproducing our results using non-social stimuli as CSs. Another limitation is the use of a twin sample, which may affect both the generalizability and validity of our findings. Regarding generalizability, small to moderate differences between twins and the normal population have been reported for some measures of fear and anxiety (Munn-Chernoff et al., 2013;Kendler et al., 1995) but not others (Pulkkinen et al., 2003).
To what extent twins diverge from the normal population with regards, specifically, to neural and physiological responding during fear conditioning is largely unknown, and, while we know of no reason to assume a significant difference between the two populations, ideally our results should be extended to a random sample from the normal population. Furthermore, regarding the validity of our findings, there is evidence demonstrating moderate heritability of SCR (Hettema et al., 2003) and BOLD-responses (Kastrati et al., 2022) during fear conditioning, meaning that twin pairs may be more similar to each other than expected by chance and therefore making data points less independent than normally assumed. To what extent this affects the reported associations to neural activity is difficult to determine. In order to test for any biases, we repeated our analysis using only single-twins, consequently cutting our sample size in half (see Appendix 7). All implicated regions evidenced correlation coefficient strengths similar, or only slightly lower, than those resulting from use of the full sample. This indicated that our results from the full sample were also applicable to unrelated individuals (see Appendix 7). In summary, while the use of a twin sample is a limitation of the present study, we contend that there is no reason to conclude that it invalidates the reported results. Finally, it should be noted that only a few studies to date have examined the longitudinal (test-retest) reliability of individual differences in discriminatory fear learning measured by SCR (Cooper et al., 2022, Preprint;Klingelhöfer-Jens et al., 2022, Preprint). While one study found evidence of fair within-person stability in 51 participants across a 9-day period (Cooper et al., 2022, Preprint), another study found evidence of poor individual-level reliability in 120 participants across a 6-month period (Klingelhöfer-Jens et al., 2022, Preprint). Similarly, previous studies examining the individual-level reliability of fMRI BOLD responding during fear conditioning have reported low to moderate reliability (Ridderbusch et al., 2021;Klingelhöfer-Jens et al., 2022, Preprint). Taken together, this means that it is currently unclear to what extent our results reflect stable individual traits, as opposed to participants' particular states at the time of measurement, which limits the interpretation of our findings. In line with Klingelhöfer-Jens et al., 2022 (Preprint) we encourage future research on individual differences in fear conditioning to explore new ways of improving the reliability of measurement.
In summary, the present study is the largest study to date on the neural correlates of autonomic fear conditioning and identified several novel areas whose activations predict individual differences in conditioned SCR. Previous findings from smaller studies could also be confirmed. Our results are consistent with the activation of a large-scale midcingulo-insular network substantiating autonomicinteroceptive processing and salience detection. We propose that altered activity within this network underlies individual differences in conditioned SCR and possibly autonomic regulation in pathological fear and anxiety. Future research should continue investigating this network as well as its possible relationship to fear and anxiety. Ultimately, such efforts may uncover the mechanisms of fear-related psychopathology and aid in its treatment.

Materials and methods Participants
This study was part of a twin study of genetic influences on fear-related brain functions. Results describing genetic influences on SCR and fMRI responses during fear conditioning will be reported elsewhere. A total of 311 adult monozygotic (N = 138) and dizygotic (N = 147) twin volunteers were recruited from the Swedish Twin Registry (Svenska Tvillingregistret). Six participants were excluded before data collection because they were unable to undergo magnetic resonance imaging. After data collection, another twenty participants were excluded from analysis due to one or several of the following reasons: unsuccessful recording of skin conductance responses (2 participants); loss of brain imaging data due to excessive head-movement (5 participants); failure to indicate shock expectancy with button press in at least 80% of trials (11 participants; see section 4.3.1); use of psychotropic medication (7 participants). Thus, 285 participants (female = 167, mean age = 33.92 years, age range = 20-58 years, SD = 10.11 years) were included in the analyses. All of these 285 participants (138 monozygotic and 147 dizygotic twins) passed the following exclusion criteria: pregnancy, inability to lie still for a 1 hr duration, intolerance of tight confinements, ongoing psychological treatment, metal objects in the body (due to surgery, fragmentation, etc.), ongoing substance abuse, or use of psychotropic medications. Non-psychotropic medication was not an exclusion criterion. However, in order to ensure the reliability of our findings, we also performed an additional supplementary analysis wherein all participants with fMRI and SCR data were included (N = 303; see Appendix 4). Participants receiving psychological treatment remained excluded as treatment could affect brain responses to emotional stimuli. Although this may be less problematic for the current analysis of brain correlates of SCR, it could be a problem when analyzing correlations between members of a twin pair if one individual receives treatment and one does not (these results will be reported elsewhere). Participants provided written informed consent in accordance with guidelines from the Regional Ethical Review Board in Uppsala and received SEK 1000 as reimbursement for their participation. The study protocol was approved by the Regional Ethical Review Board in Uppsala (Dnr 2016/171).

Stimuli and contexts
During the discriminative fear conditioning task, two male, three-dimensional, virtual humanoid characters and a virtual environment (Figure 3) were used. The characters and the environments were created in Unity (version 5.2.3, Unity Technologies, San Francisco, CA) and consisted of a room with four red brick walls, a grey concrete roof, and a wooden floor.

Stimulus presentation software
The virtual characters and the environment were shown on a flat surface in the MR scanner by a projector (Epson EX5260). Stimulus presentation was handled by the Unity 3D-engine (version 5.2.3, Unity Technologies, San Francisco, CA). The Unity software communicated with BIOPAC (BIOPAC Systems, Goleta, CA) through a custom-made serial interface using standard libraries by Microsoft (Microsoft Corporation, Albuquerque, New Mexico).

Brain imaging
Imaging data were acquired using a 3.0 T scanner (Discovery MR750, GE Healthcare) and an eightchannel head-coil. Foam wedges, earplugs, and headphones were used to reduce head motion and scanner noise. We acquired T1-weighted structural images with whole-head coverage, TR = 6400ms, TE = 28ms, acquisition time = 6.04 min, and flip angle = 11°. Functional images were acquired using gradient echo-planar-imaging (EPI), TR = 2400ms, TE = 28ms, flip angle = 80°, slice thickness = 3.0 mm with no spacing, axial orientation, frequency direction R/L, interleaved (bottom up), number of slices were 47 and voxel size 3.0 mm 3 . Higher order shimming was performed, and five dummy scans were run before the experiment.

Skin conductance responses
Skin conductance was recorded with the MP-150 BIOPAC system (BIOPAC Systems, Goleta, CA). Radio-translucent disposable dry electrodes (EL509, BIOPAC Systems, Goleta, CA) were coated with isotonic gel (GEL101, BIOPAC Systems, Goleta, CA) and placed on the palmar surface of the participant's left hand. The signal was high-pass (0.05 Hz) filtered using the built-in BIOPAC hardware Butterworth filter. SCRs were scored using Ledalab software package (v 3.4.9) (Benedek and Kaernbach, 2010) implemented in Matlab 2020a (Mathworks, Inc, Natick, MA). Minimum response threshold was set to 0.01 µS. After filtering and before analysis, the SCR signal was down-sampled from 2000 Hz to 200 Hz (factor mean). The options specified for the Ledalab batch run were 'open', 'biotrace', 'downsample', 10, 'analyze','CDA', 'optimize',4. SCR was analyzed using standard peak score (throughto-peak, TTP.AmpSum) 1-4 s after CS onset for each participant. To check that whole brain SCR correlations were not dependent on the choice of peak scoring window, we also analyzed SCR with a window of 1-5 s after CS onset. We also scored SCR using a software package called PsPM (Bach and Friston, 2013), which uses a model-based approach in estimating SCR (see Appendix 8 for details). We performed these variants of SCR scoring as part of a sensitivity analysis to ensure that correlation results between SCR and brain activity were not dependent on the choice of SCR scoring method.
Unlike previous studies considering the neural correlates of CS+ > CS- SCR (LaBar et al., 1998;Phelps et al., 2004;Petrovic et al., 2008;MacNamara et al., 2015;Marin et al., 2020), SCRs in the present study were range-corrected by Z transformation within individuals (BenShakhar, 1985). Z transformation increases sensitivity to conditioning-related effects and prevents confounding by non-conditioning-related individual differences in general SCR magnitude (BenShakhar, 1985;Staib et al., 2015). For comparison, however, correlations based on raw SCR scores were also examined and can be found in Appendix 3.
Skin conductance was recorded without a low-pass filter. By using this recording procedure, we noticed that a small number of trials produced unreasonably high SCR values (e.g. 17 mS responding), likely due to electrode movement. As such extreme values may skew correlations even using standard scores, we excluded from all analysis trials with a raw value SCR score > 3 mS. This was based on previous research indicating a general maximum SCR between 2 and 3 mS using similar methodology as the one used in this paper (Boucsein, 2012). Using this criterion 97/12200 trials were excluded from analysis (0.795% of all trials).

Fear conditioning task
Two virtual characters served as CSs, one as a threat cue (CS+) predicting the US and the other as a control cue (CS-). CSs were presented on a screen in the MR scanner at a distance of 2.7 m in the virtual environment. The relatively long distance of 2.7 m was selected in order for the effect of conditioning on SCR not to be occluded by proximal threat effects on SCR, as was observed in two previous studies by Rosén et al., 2017;Rosén et al., 2019. Participants were told prior to the experiment that they could learn to predict the US but were not told which character served as the CS+. Participants were furthermore instructed to indicate which character would be followed by the US by selecting either 'yes' or 'no' by pressing a button immediately following each CS presentation. The inclusion of the button presses was to assure participants' attention during the task as well as to confirm that participants, overall, learned the contingency. Which of the two characters that served as CS+ and CS-was counterbalanced across participants using four different stimulus presentation orders. Each CS presentation lasted for 6 s followed by an inter-stimulus interval of 8-12 s, during which the context was still displayed, but no CS was present.
A habituation phase preceded the fear conditioning, in which each CS was presented four times without reinforcement for a total of eight CS presentations. This was followed by the fear conditioning phase, during which each CS type (CS+ and CS-) was presented 16 times. The experimental task thus consisted of 40 CS presentations in total: 8 during habituation and 32 during conditioning. During conditioning, eight of the CS+ presentations co-terminated with a presentation of the US (50% partial reinforcement schedule; in accordance with guidelines for increasing sensitivity to inter-individual differences in fear conditioning, see Lonsdorf and Merz, 2017). Total time for the fear conditioning task was 9 min and 47 s.
The US consisted of an electric shock delivered to the subjects' wrist via radio-translucent dry electrodes (EL509, BIOPAC Systems, Goleta, CA). Prior to the experiment, the shock was calibrated using an ascending staircase procedure wherein shock intensity is increased until rated by participants as 'uncomfortable' but not 'painful' (Åhs et al., 2015). US duration was 16ms and controlled using the STM100C module connected to the STM200 constant voltage stimulator (BIOPAC Systems, Goleta, CA).
In all sequences, the first CS+ presentation following the 4 CS+ habituation trials was always reinforced. The sequences differed in whether the CS-or CS+ started the acquisition phase. If the reinforced CS+ is always the first trial in the acquisition phase, the CS-trial following the US will be elevated due to sensitization. This was why the presentation order was counterbalanced.
Analysis of SCR data SCR Z scores were averaged separately across CS+ and CS-trials within each participant. A paired samples t test was performed to compare the average CS+ SCRs to the average CS-SCRs at an alpha level (of significance) of α = .05 using JASP software (version 0.14.1, JASP Team (2020)). This allowed us to determine whether the fear conditioning task was successful in evoking greater SCR to the CS+ than to the CS-. Secondly, in order to examine the correlations between conditioned SCR and fMRI responses during fear conditioning, an SCR difference score was calculated for each participant by subtracting the average SCR to CS+ presentations from the average SCR to CS-presentations. The distribution of SCR difference scores was examined to ensure the validity and sensitivity of neural regression analyses (see Appendix 1 regarding methodology and results of SCR distribution analysis).

Online recording of shock expectancy
Participants pressed one of two buttons each time a CS was displayed to indicate whether they were expecting to receive an electric shock (coded 1) or not (coded 0). The mean response was computed for the CS+ and the CS-presentations. A t-test was performed in JASP to compare mean shock expectancy to the CS+ and CS-, indicating whether participants learned the contingency.

Analysis of fMRI data
Analyses of fMRI data were performed using SPM12 (version 6685, Wellcome Department of Cognitive Neurology, University College, London). Preprocessing of images included interleaved slice time correction and realignment of functional volumes. For each participant, the mean functional image was co-registered to the anatomical T1-weigthed image. Quality control of functional images was performed using MRIQC v0.16 (mask validation matching; Esteban et al., 2017). Realignment parameters were inspected for excessive movement (defined as 5 mm) during scanning. Anatomical images were segmented using 4 tissue classes and normalized to Montreal Neurological Institute (MNI) standard space. The co-registered functional images were next warped to MNI space using the same parameters that were used for the anatomical image. An 8 mm FHMW Gaussian kernel was used for smoothing of the functional images.
First-level analysis used event-related modeling including regressors for CS+ habituation trials, CS-habituation trials, CS+ acquisition trials, CS-acquisition trials, and US delivery in a general linear model (GLM). Regressors mapped to the intervals of 6 s where each type of stimulus was displayed in the scanner and was convolved with the hemodynamic response function to predict the fMRI time course (for the brief US, the regressor still mapped to the 6 s following US delivery). Also, 6 movement parameter regressors, derived from the image realignment, and a mean value intercept regressor (a vector of ones) were included in the model. At the second level, we first examined the overall whole brain CS+ > CS-BOLD contrast to assess conditioning-related effects. This analysis included both reinforced and non-reinforced CS+ trials, as the US was delivered 2 s after the 4 s sampling window for SCR scoring, and therefore could not confound the SCRs on reinforced trials. Voxel-wise statistical significance was calculated using t tests implemented in the SPM12 software with an alpha level of α = .05 using family-wise error (FWE) correction. As results regarding this contrast are not central to the present study, they are published in Appendix 6. However, it should be noted that these results were typical for fear conditioning studies in general (see e.g. Fullana et al., 2016). In the present analysis, we proceeded to examine regional activity within the whole brain CS+ > CS-BOLD contrast that covaried with individual differences in the SCR difference score. This was done by entering each participant's previously obtained differential SCR score (CS+ minus CS-) as a second level betweensubjects regressor of their average CS+ > CS-BOLD activation, effectively correlating the CS+ > CS-SCR difference with the CS+ > CS-BOLD difference. Voxel-wise statistical significance was again calculated using t-tests implemented in the SPM12 software with an FWE corrected alpha level of α = .05. Notably, this provided a local maxima height threshold of t = 4.36 given our sample size (N = 285). This roughly corresponds to a correlation coefficient of r = .25 given the conversion formula r = sqrt (t 2 /(t 2 + DF)), thus corresponding to an effect size within the medium range according to guidelines by Cohen, 1988;Cohen, 1992. Such an effect size is lower than estimated effect sizes from the study by MacNamara et al., 2015, the only previous study meeting the minimum requirements for sample size suggested by Yarkoni and Braver, 2010 and performing a correlation analysis similar to the present one. Specifically, by using the conversion formula r = sqrt (z 2 /(z 2 +N)), we estimated effect sizes in this study to be large (r = .48 and r = .49). Based on this observation, the power of the present study was considered acceptable.
Secondly, we tested a hypothesized association between individual differences in amygdala response and SCR using a region-of-interest (ROI) analysis. The amygdala ROI was defined using the automated anatomical labeling (AAL version 1) library (Tzourio-Mazoyer et al., 2002) and included both right and left amygdala. Analysis was performed the same way as the previous whole brain correlation except that individual SCR scores were this time correlated exclusively to contrast values within the amygdala ROI, thus increasing sensitivity within this theoretically implicated region. This analysis yielded a local maxima height threshold of t = 3.66 given our sample size (N = 285). This roughly corresponds to a correlation coefficient of r = 0.21 given the conversion formula r = sqrt (t 2 / (t 2 +DF)), thus corresponding to an effect size within the small to medium range according to guidelines by Cohen, 1988;Cohen, 1992. Third, in order to compare the potentially independent contributions to individual differences in conditioned SCR of different clusters of voxels that showed significant correlations in the whole-brain analysis, we extracted the eigenvariates of the first-level CS+ > CS-contrast values in each cluster of voxels. The SPM12 software's built-in 'extract eigenvariates' function was used to extract eigenvariates from unadjusted contrast values. SPM12 uses Singular Value Decomposition (SVD) to calculate eigenvariates. The result was one value per cluster of voxels (see Table 2 for the list of regions) for each individual. While eigenvariates are strongly correlated (r > .9) to the mean contrast value of the same voxels, they are used instead of the mean because they are less sensitive to extreme values in individual voxels (for more thorough explanations, see Ridgeway, 2012;Penny, 2004). The vectors of eigenvariates for each region were next used as a regressors of differential SCR scores in a hierarchical regression analysis implemented in the JASP software. The potentially unique contribution to SCR of each region could then be tested by examining the significance of the beta weights of each region within the model.
Finally, in order to determine if the obtained neural correlates to SCR difference scores were driven more by SCR to CS+ or CS-, we correlated eigenvariates extracted from implicated neural regions and average SCRs to the CS+ and CS-separately. In this analysis, we used square root transformed raw value SCRs instead of Z-transformed SCR in order to obtain roughly normalized SCRs without confounding of CS+ and CS-response magnitude, such as that which occurs when using Z transformation. Specifically, since the Z transformed SCR is defined as the difference between the SCR on a given trial (or trial type) and the average SCR across all trials, divided by the standard deviation of SCRs across all trials, this means that using Z scores conflates CS+ and CS-responding with responses to both average CS+ and CS-inclusively, and therefore cannot be used to determine the influence of neural activity on a specific trial type. As the distributions of square root transformed raw value SCRs to the CS+ and CS-still did not meet criteria for normality (by visual inspection of histograms, QQ-plot, and Shapiro-Wilke's test showing p < 0.001; see Appendix 2), we used Spearman's Rho instead of Pearson's r correlations. To compare the difference between CS+ and CS-correlations we used the Steiger, 1980 direct comparison of dependent correlation coefficients as implemented in free automated software by Lee and Preacher, 2013, as this is the most robust way of testing the difference between dependent Spearman coefficients (Myers and Sirois, 2006). To compensate for multiple testing during correlation comparisons, we used a Bonferroni corrected alpha level of α = 0.05/33 = 0.00151. This compensated for a total of 33 tests, reflecting 11 implicated neural regions each being correlated to CS+ and CS-separately, as well as each having these correlations compared in an additional test. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Additional files
Supplementary files • Transparent reporting form

Data availability
The present study reports data from human participants that did not explicitly consent to their raw neuroimaging and physiological data being made public. Therefore, raw neuroimaging and physiological data from the present study cannot currently be made publicly available. Requests for the anonymized raw neuroimaging and physiological data should be made to Fredrik Åhs ( fredrik. ahs@ miun. se) and will be reviewed by an independent data access committee, taking into account the research proposal and the intended use of the data. Requestors are required to sign a data transfer agreement to ensure participants' confidentiality is maintained prior to release of any data and that procedures conform with the EU legislation on the general data protection regulation and local ethical regulations. While access to raw source data is thus limited, processed data and standardized statistical images sufficient to reproduce the reported results and figures are publicly and freely available at https://osf.io/7dz9p/. Specifically, we provide statistical brain images in NIfTI file format used to render Figure 1a, Figure 2a, Table 2, Appendix 3-figure 1, Appendix 4- figure 1 and Appendix 6figure 1 of the present study. We also provide brief explanations of the software used to produce all source data files, along with the SPM job files used for neuroimaging analyses. In the event that ethical approval to publicly share the raw neuroimaging data of the present study is obtained at a later stage; this data will also be made publicly available on the OSF site. In the present journal we have included Figure 1-source data 1, which provides source data for Figure 1b, Figure 2b and statistical analyses reported in section 2.2.3 as well as for Appendix 1-figure 1, Appendix 2-figure 1 and Appendix 5-table 1.
The following dataset was generated: Note. MNI coordinates and t values represent significant peak voxels of each cluster. Statistical significance was calculated using t tests implemented within the SPM software with an FWE corrected alpha level of α = .05.
In summary, both supplementary whole brain analyses (Appendix 3 and 4) implicated largely the same set of regions as the main analysis of the main text (see Table 2). In particular, rightlateralized regional activations in the dorsal anterior cingulate cortex/anterior midcingulate cortex, anterior insula, inferior frontal gyrus, temporoparietal junction, superior frontal gyrus/ dorsal premotor cortex and midbrain were consistent across analyses. Thus, findings regarding these regions appear robust. However, neural activations in the left superior frontal gyrus/dorsal premotor cortex, left temporoparietal junction, left inferior frontal gyrus and right superior parietal lobe were not consistently implicated. Hence, findings regarding these latter regions do not appear as robust. However, two important things should be noticed regarding these latter regions. First, our conclusions regarding the main findings of the study are not heavily dependent on the activation of these latter regions. Second, for reasons explained in the Materials and Methods section of the main text, we consider the main analysis of the main text to be the most valid analysis. intolerance of tight confinements, ongoing psychological treatment, metal objects in the body (due to surgery, fragmentation etc.), ongoing substance abuse, use of psychotropic medications, unsuccessful recording of skin conductance responses, loss of brain imaging data due to excessive head movement, participant failure to comply with task instruction regarding button press in at least 80% of trials. Note. MNI coordinates and t values represent significant peak voxels within each cluster. Statistical significance was calculated using t tests with an FWE corrected alpha level of α = .05 within the SPM software.
time window following CS-onset for peak-detection. Non-reinforced trials only included non-reinforced CS+ trials when computing SCR and neural responses. To control for shock expectancy effects, we included average shock expectancy (rated online as 0 or 1) as a covariate correlating with whole-brain responses. We also controlled for eventual familial influences by splitting twin pairs and correlating SCR in the two samples (first twin, second twin) to brain contrast values. Note that the full sample included twin pairs as well as twins without a sibling, hence there is a discrepancy in the number of twins for first and second columns (135 vs 146).
Correlations between SCR to US and CS+ > CS-contrast-values are shown in the last column for reference.

Anatomical region
Hemisphere Voxels Note. MNI coordinates and t values represent significant peak voxels of clusters from the main analysis. Statistical significance was calculated using t tests implemented within the SPM software with an FWE corrected alpha level of α = .05.