## Abstract

During real-life situations, multiple factors interact dynamically to determine threat level. In the current fMRI study involving healthy adult human volunteers, we investigated interactions between proximity, direction (approach vs. retreat), and speed during a dynamic threat-of-shock paradigm. As a measure of threat-evoked physiological arousal, skin conductance responses were recorded during fMRI scanning. Some brain regions tracked individual threat-related factors, and others were also sensitive to combinations of these variables. In particular, signals in the anterior insula tracked the interaction between proximity and direction where approach versus retreat responses were stronger when threat was closer compared with farther. A parallel proximity-by-direction interaction was also observed in physiological skin conductance responses. In the right amygdala, we observed a proximity by direction interaction, but intriguingly in the opposite direction as the anterior insula; retreat versus approach responses were stronger when threat was closer compared with farther. In the right bed nucleus of the stria terminalis, we observed an effect of threat proximity, whereas in the right periaqueductal gray/midbrain we observed an effect of threat direction and a proximity by direction by speed interaction (the latter was detected in exploratory analyses but not in a voxelwise fashion). Together, our study refines our understanding of the brain mechanisms involved during aversive anticipation in the human brain. Importantly, it emphasizes that threat processing should be understood in a manner that is both context-sensitive and dynamic.

## INTRODUCTION

Anticipation of aversive events leads to a repertoire of changes in behavioral, physiological, and brain responses that contribute to the handling of the negative consequences of such events. At the same time, abnormalities in aversive anticipatory processing are thought to underlie many mental disorders, such as anxiety and depression (Dillon et al., 2014; Grupe & Nitschke, 2013). Hence, understanding the brain mechanisms of aversive anticipation is important from both basic and clinical standpoints.

In humans, aversive anticipation has been investigated with paradigms in which punctate cues signal an upcoming negative event (Brown, Seymour, Boyle, El-Deredy, & Jones, 2008; Nitschke, Sarinopoulos, Mackiewicz, Schaefer, & Davidson, 2006; Simmons, Strigo, Matthews, Paulus, & Stein, 2006; Jensen et al., 2003; Böcker, Baas, Kenemans, & Verbaten, 2001) or by blocked manipulations with constant threat level (McMenamin, Langeslag, Sirbu, Padmala, & Pessoa, 2014; Vytal, Overstreet, Charney, Robinson, & Grillon, 2014). However, during most real-world situations, aversive anticipation changes dynamically over time. An important factor in determining threat level is proximity, as when a prey reacts differently to the presence of a predator when the latter is proximal compared with distant (Figure 1A; Blanchard, Griebel, Pobbe, & Blanchard, 2011; Blanchard & Blanchard, 1990). Other factors involve direction, namely whether threat is approaching versus retreating (Figure 1B) and speed, reflecting how fast or slow the threat is moving (Fanselow & Lester, 1988). Some studies have taken initial strides at investigating how some of these factors influence brain responses during aversive anticipation. For instance, the contrast of proximal versus distal threats revealed fMRI responses in a host of brain regions, including the anterior insula, midbrain periaqueductal gray (PAG), and bed nucleus of the stria terminalis (BST; Mobbs et al., 2010; Somerville, Whalen, & Kelley, 2010); evidence for amygdala involvement linked to threat proximity is mixed (Mobbs et al., 2010; Somerville et al., 2010). Similarly, comparison of approaching versus retreating threats has revealed responses in the anterior insula, BST, and amygdala (Mobbs et al., 2010).

Figure 1.

Threat-related factors and their interaction. (A) Closer and farther threat, where threat is represented by an aversive shock when circles touched. (B) Direction of threat: approach versus retreat. (C) Threat level may depend on both proximity (closer and farther) and direction (left panels indicate approach; right panels indicate retreat).

Figure 1.

Threat-related factors and their interaction. (A) Closer and farther threat, where threat is represented by an aversive shock when circles touched. (B) Direction of threat: approach versus retreat. (C) Threat level may depend on both proximity (closer and farther) and direction (left panels indicate approach; right panels indicate retreat).

Thus far, studies have considered the effects of threat proximity and direction independently. Hence, it is currently unknown how such factors potentially interact in the brain during aversive anticipation (Figure 1C). This is an important gap in our knowledge base because behavioral findings have extensively documented interactions between threat-related factors, which have produced several influential theoretical accounts (for excellent discussion, see Mobbs, Hagan, Dalgleish, Silston, & Prévost, 2015). Furthermore, it is not only important to investigate how multiple threat-related factors interact but to understand how the brain tracks them continuously. In particular, do signal fluctuations in brain regions track threat-related factors dynamically? If so, to what factor(s) and factor combinations are they sensitive?

To address these questions, we devised a paradigm in which threat was dynamically modulated during fMRI scanning. Two circles moved on the screen, sometimes moving closer and sometimes moving apart, and at varying speeds (Figure 2). Participants were instructed to pay attention to the circles on the screen and were explicitly informed that, if they touched, the participants would receive an unpleasant shock. As a measure of threat-evoked physiological arousal, skin conductance responses (SCRs) were recorded during scanning. Our paradigm allowed us to investigate the role played by the interaction between proximity (nearer vs. farther circles), direction (approach vs. retreat), and speed (faster vs. slower) in determining brain responses during anticipatory threat processing. Importantly, the impact of the factors “proximity” and “speed” were assessed parametrically (i.e., continuously) as they varied dynamically. Therefore, the paradigm allowed us to test how multiple threat-related factors “dynamically” influence signals fluctuations across brain regions. Specifically, do they provide independent contributions or do they interact in regions important for threat processing, such as the anterior insula, amygdala, PAG, and BST? Intuitively, probing interactions allowed us to evaluate the extent to which the influence of one factor on threat anticipation depended on the values of other factor(s). For instance, in terms of a two-way interaction, we anticipated that the influence of direction (i.e., approaching vs. retreating threat) would depend on proximity (i.e., whether the threat was near versus far; Figure 1C). In terms of three-way interactions, we sought to evaluate if the interaction between the continuously manipulated factors of proximity and speed depended on direction.

Figure 2.

Experimental paradigm. Two circles moved randomly on the screen and a shock was administered to the participant if they touched. The inset represents threat proximity (the distance between the two circles), which varied continuously. A central goal of the study was to determine the extent to which signal fluctuations in brain regions (such as the anterior insula) followed threat-related factors (including proximity) and their interactions.

Figure 2.

Experimental paradigm. Two circles moved randomly on the screen and a shock was administered to the participant if they touched. The inset represents threat proximity (the distance between the two circles), which varied continuously. A central goal of the study was to determine the extent to which signal fluctuations in brain regions (such as the anterior insula) followed threat-related factors (including proximity) and their interactions.

## METHODS

### Participants

Eighty-five participants (41 women, ages 18–40 years, average = 22.62 years, SD = 4.85) with normal or corrected-to-normal vision and no reported neurological or psychiatric disease were recruited from the University of Maryland community (of the original sample of 93, data from seven participants were discarded because of technical issues during data transfer [specifically, field maps were lost], and one other participant was removed because of poor structural–functional alignment). The project was approved by the University of Maryland, College Park institutional review board, and all participants provided written informed consent before participation. The data analyzed here were investigated in an entirely separate fashion at the level of networks and published previously (Najafi, Kinnison, & Pessoa, 2017). The sample size was not based on an explicit statistical power analysis. At the outset, we sought to collect around 90 participants to allow investigation of the data in terms of separate “exploratory” and “test” sets in the network study (Najafi et al., 2017). For the investigation of activation (this article), our intention was to use the available data in a single type of analysis.

### Anxiety Questionnaires

Participants completed the trait portion of the Spielberger State–Trait Anxiety Inventory (Spielberger, Gorsuch, & Lushene, 1970) before scanning (average = 17.23 days, SD = 15.90) and then completed the state portion of the State–Trait Anxiety Inventory immediately before the scanning session.

### Procedure and Stimuli

Two circles with different colors moved around on the screen randomly. When they collided with each other, an unpleasant mild electric shock was delivered. Overall, proximity, direction of movement, and relative speed of the circles were used to influence perceived threat. The position of each circle (on the plane), xt, was defined based on its previous position, xt−1, plus a random displacement, Δxt:
$xt=xt−1+Δxt$
The magnitude and direction of the displacement was calculated by combining a normal random distribution with a momentum term to ensure motion smoothness, while at the same time remaining (relatively) unpredictable to the participants. Specifically, the displacement was updated every 50 msec as follows:
$Δxt=1−cΔxt−1+cN01$
where c = 0.2 and N(0, 1) indicates the normal distribution with zero mean and standard deviation of 1.

Visual stimuli were presented using PsychoPy (www.psychopy.org/) and viewed on a projection screen via a mirror mounted to the scanner's head coil. Each participant viewed the same sequence of circle movements. The total experiment included six runs (457 sec each), each of which had six blocks (3 of 85 participants had only five runs). In each block, the circles appeared on the screen and moved around for 60 sec; blocks were separated by a 15-sec off period during which the screen remained blank. Each run ended with a 7-sec blank screen.

To ensure that the effects of threat proximity and direction were uncorrelated, half of the blocks in each run were temporally reversed versions of the other blocks in that run. Temporally reversing the stimulus trajectories guarantees that proximity and direction are uncorrelated because reversing time changes the sign of the direction (i.e., approach becomes retreat). To optimize the experimental design, 10,000 candidate stimuli trajectories and block orders were generated. We then selected six runs, which minimized collinearity between all predictors of interest (see below), measured as the sum of respective variance inflation factors (Neter, Kutner, Nachtsheim, & Wasserman, 1996).

In each run, the circles collided eight times within four of six blocks (one to three times in a block); in the remaining two blocks, there were no collisions. Each collision resulted in the delivery of an electric shock. The 500-msec electric shock (composed of a series of current pulses at 50 Hz) was delivered by an electric stimulator (Model E13-22 from Coulbourn Instruments, Whitehall, PA) to the fourth and fifth fingers of the nondominant left hand via MRI-compatible electrodes. To calibrate the intensity of the shock, each participant was asked to choose his or her own stimulation level immediately before functional imaging, such that the stimulus would be “highly unpleasant but not painful.” After each run, participants were asked about the unpleasantness of the stimulus to recalibrate shock strength, if needed. SCR data were collected using the MP-150 system (BIOPAC Systems, Inc., Goleta, CA) at a sampling rate of 250 Hz by using MRI-compatible electrodes attached to the index and middle fingers of the nondominant left hand. Because of technical problems and/or experimenter errors during data collection, SCR data were not available in two participants, and six participants had only five runs of the SCR data; one participant who had only three runs of data was excluded from the analysis of SCR data.

### MRI Data Acquisition

Functional and structural MRI data were acquired using a 3-T Siemens TRIO scanner with a 32-channel head coil. First, a high-resolution T2-weighted anatomical scan using Siemens's SPACE sequence (0.8 mm isotropic) was collected. Subsequently, we collected 457 functional EPI volumes in each run using a multiband scanning sequence (Feinberg et al., 2010), with repetition time = 1.0 sec, echo time = 39 msec, field of view = 210 mm, and multiband factor = 6. Each volume contained 66 nonoverlapping oblique slices oriented 30° clockwise relative to the AC–PC axis (2.2 mm isotropic). A high-resolution T1-weighted MPRAGE anatomical scan (0.8 mm isotropic) was collected. In addition, in each session, double-echo field maps (TE1 = 4.92 msec, TE2 = 7.38 msec) were acquired with acquisition parameters matched to the functional data.

### fMRI Preprocessing

To preprocess the fMRI and anatomical MRI data, a combination of packages and in-house scripts were used. The first three volumes of each functional run were discarded to account for equilibration effects. Slice-timing correction—with Analysis of Functional Neuroimages' (AFNI; Cox, 1996) 3dTshift—used Fourier interpolation to align the onset times of every slice in a volume to the first acquisition slice, and then a six-parameter rigid body transformation (with AFNI's 3dvolreg) corrected head motion within and between runs by spatially registering each volume to the first volume.

In this study, we strived to improve functional–anatomical coregistration given the small size of some of the structures of interest. Skull stripping determines which voxels are to be considered part of the brain and, although conceptually simple, plays a very important role in successful subsequent coregistration and normalization steps. Currently, available packages perform suboptimally in specific cases, and mistakes in the brain-to-skull segmentation can be easily identified. Accordingly, to skull strip the T1 high-resolution anatomical image (which was rotated to match the oblique plane of the functional data with AFNI's 3dWarp), we used six different packages, ANTs (Avants, Tustison, & Song, 2009; http://stnava.github.io/ANTs/), AFNI (Cox, 1996; http://afni.nimh.nih.gov/), ROBEX (Iglesias, Liu, Thompson, & Tu, 2011; https://www.nitrc.org/projects/robex), FSL (Smith et al., 2004; http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/), SPM (www.fil.ion.ucl.ac.uk/spm/), and BrainSuite (Shattuck & Leahy, 2002; http://brainsuite.org/), and used a “voting scheme” as follows: Based on T1 data, a voxel was considered to be part of the brain if four of six packages estimated it to be a brain voxel; otherwise, the voxel was not considered to be brain tissue (for six participants whose T1 data were lost due to issues during data transfer, the T2 image was used instead and only the ANTs package was used for skull stripping).

Subsequently, FSL was used to process field map images and create a phase distortion map for each participant (by using bet and fsl_prepare_fieldmap). FSL's epi_reg was then used to apply boundary-based coregistration to align the unwarped mean volume registered EPI image with the skull-stripped anatomical image (T1 or T2), along with simultaneous EPI distortion correction (Greve & Fischl, 2009).

Next, ANTS was used to estimate a nonlinear transformation that mapped the skull-stripped anatomical image (T1 or T2) to the skull-stripped MNI152 template (interpolated to 1-mm isotropic voxels). Finally, ANTS combined the nonlinear transformations from coregistration/unwarping (from mapping mean functional EPI image to the anatomical T1 or T2) and normalization (from mapping T1 or T2 to the MNI template) into a single transformation that was applied to map volume-registered functional volumes to standard space (interpolated to 2-mm isotropic voxels). In this process, ANTS also utilized the field maps to simultaneously minimize EPI distortion. The resulting spatially normalized functional data were blurred using a 4-mm FWHM Gaussian filter. Spatial smoothing was restricted to gray matter mask voxels. Finally, the intensity of each voxel was normalized to a mean of 100 (separately for each run).

### Voxelwise Analysis

Each participant's preprocessed fMRI data were analyzed using multiple linear regression with AFNI (restricted to gray matter voxels) using the 3dDeconvolve program (https://afni.nimh.nih.gov/afni/doc/manual/3dDeconvolve.pdf). Time series data were analyzed according to the following model (additional nuisance variables are described below):
$Y=βPP+βDD+βSS+βPDPD+βPSPS+βDSDS+βPDSPDS$
(1)
where P indicates proximity, D represents direction, and S represents speed. Variables were determined based on circle positions on the screen. Proximity was defined as the Euclidean distance between the two circles; direction indicated approach versus retreat; speed was the discrete temporal difference of proximity. The products PD, PS, and PDS represent the interactions terms; the individual terms P, D, and S were mean-centered before multiplication to reduce potential collinearity. The resulting regressors exhibited pairwise correlations that were relatively small (the largest was .41), and all variance inflation factors were less than 1.3, indicating that model estimation was unproblematic (Mumford, Poline, & Poldrack, 2015).

In addition to the variables above, we included regressors for visual motion (velocity tangential to the difference vector of the combined circle-to-circle stimulus), sustained block event (60-sec duration), and block-onset and block-offset events (1-sec duration) to account for transient responses at block onset/offset. All regressors were convolved with a standard hemodynamic response based on the gamma variate model (Cohen, 1997). Note that interaction regressors were multiplied before convolution; also, as stimulus-related display information was updated every 50 msec (20 Hz), convolution with the hemodynamic response was performed before decimating the convolved signal to the fMRI sample rate (1 Hz). To simplify plotting, decimated regressors were scaled by their corresponding root mean square value (thus, multiplicative interactions terms were on the same scale as simple effects). Other regressors included in the model included six motion parameters (three linear displacements and three angular rotations) and their discrete temporal derivatives. To further control for head motion-related artifacts in the data (Siegel et al., 2014), we excluded volumes (on average 0.4%) with a frame-to-frame displacement of more than 1 mm. To model baseline and drifts of the MRI signal, regressors corresponding to polynomial terms up to the fourth order were included (for each run separately). Finally, to minimize effects due to the physical shock event, data points in a 15-sec window after shock delivery were discarded from the analysis. It should be pointed out that to partly account for the fact that the circles were most proximal just before shock events, the design included time periods when circles were very close but did not touch eventually.

### Group Analysis

Whole-brain voxelwise random-effects analyses were conducted using response estimates from individual-level analyses (restricted to gray matter voxels) in AFNI. To probe the effects of the regressors of interest, we ran separate one-sample t tests against zero using the AFNI's 3dttest++ program.

The alpha-level for voxelwise statistical analysis was determined by simulations using the 3dClustSim program (restricted to gray matter voxels). For these simulations, the smoothness of the data was estimated using 3dFWHMx program (restricted to gray matter voxels) based on the residual time series from the individual-level voxelwise analysis. Taking into account the recent report of increased false-positive rates linked to the assumption of Gaussian spatial autocorrelation in fMRI data (Eklund, Nichols, & Knutsson, 2016), we used the -acf (i.e., autocorrelation function) option recently added to the 3dFWHMx and 3dClustSim tools, which models spatial fMRI noise as a mixture of Gaussian plus monoexponential distributions. This improvement was shown to control false-positive rates around the desired alpha level, especially with relatively stringent voxel-level uncorrected p values such as .001 (Cox, Chen, Glen, Reynolds, & Taylor, 2017). Based on a voxel-level uncorrected p value of .001, simulations indicated a minimum cluster extent of 13 voxels (2.0 × 2.0 × 2.0 mm) for a cluster-level corrected alpha of .05.

### BST ROI Analysis

The BST is a basal forebrain region and has been frequently implicated in threat-related processing (Fox, Oler, Tromp, Fudge, & Kalin, 2015; Davis, Walker, Miles, & Grillon, 2010), along with other regions such as the amygdala and anterior insula (Pessoa, 2016). Because the BST is a small region, analysis based on spatially smoothed data would be susceptible to signals from surrounding structures. To reduce this possibility, we conducted an additional BST ROI analysis using spatially unsmoothed data. Bilateral BST ROIs were defined anatomically according to the probabilistic mask of the BST (at 25% threshold) recently reported by Blackford and colleagues (Theiss, Ridgewell, McHugo, Heckers, & Blackford, 2017). For this analysis, no spatial smoothing was applied. In each participant, for each ROI, a representative time series was created by averaging the unsmoothed time series from all the gray matter voxels within the anatomically defined ROI (left: nine voxels, right: eight voxels). Then, as in the individual-level voxelwise analysis, multiple linear regression analysis was run using the 3dDeconvolve program to estimate condition-specific responses. At the group level, as in the voxelwise analysis, we ran separate one-sample t tests against zero using the corresponding regression coefficients from the individual-level analysis.

### SCR Analysis

Each participant's SCR data were initially smoothed with a median filter over 50 samples (200 msec) to reduce scanner-induced noise. In each run, the first 3 sec of data were discarded (corresponding to first three volumes excluded in the fMRI analysis), and the remaining data were resampled by decimating the 250-Hz sample rate to the sample rate of fMRI data (1 Hz) and subsequently Z scored. The preprocessed SCR data were then analyzed using multiple linear regression using the 3dDeconvolve program in AFNI (for related approaches, see Bach, Flandin, Friston, & Dolan, 2009; Engelmann, Meyer, Fehr, & Ruff, 2015). We used the same regression model as the one used for fMRI data (see Equation 1). In addition, we included regressors for visual motion (velocity tangential to the difference vector of the combined circle-to-circle stimulus), sustained block event (60-sec duration), and block-onset and block-offset events (1-sec duration) to account for transient responses at block onset/offset. All regressors were convolved with a canonical SCR model based on the sigmoid exponential function (Lim et al., 1997; Figure 3). In addition, constant and linear terms were included (for each run separately) to model baseline and drifts of the SCR. To minimize effects due to the physical shock event, data points in a 15-sec window after shock delivery were discarded from the analysis. At the group level, to probe the effects of the regressors of interest, we ran separate one-sample t tests against zero using the corresponding regression coefficients from the individual-level analysis.

Figure 3.

Skin conductance response model based on the sigmoid exponential function (Lim et al., 1997). A.U. = arbitrary units.

Figure 3.

Skin conductance response model based on the sigmoid exponential function (Lim et al., 1997). A.U. = arbitrary units.

### Relationship between SCR and Brain Activity

To probe the relationship between brain activity and physiological arousal, we focused on the right anterior insula and the right amygdala clusters that exhibited a proximity by direction interaction (see Results). For each cluster, an interaction index was created by averaging the corresponding regression coefficients (βPD in Equation 1) from all the voxels within the cluster (after cluster-level thresholding). Then, for each cluster, we ran a robust correlation (Rousselet & Pernet, 2012; Wilcox, 2012) across participants. For each participant, we considered the average fMRI interaction regression coefficient and the corresponding interaction term in the SCR data (specifically, the coefficient βPD obtained from the SCR regression analysis).

### Relationship between Threat Anticipation and Physical Shock Responses

In an exploratory analysis, we probed the relationship between activity related to threat anticipation and responses to physical shock itself. For the anticipatory activity, we considered the proximity by direction interaction and focused on the right anterior insula and right amygdala clusters, which exhibited this interaction (see Results). To estimate responses to physical shocks, we ran a separate multiple regression analysis with all the regressors as in the original model along with an additional regressor that modeled physical shock events (500 msec). As noted above, these events were discarded in the main analyses to minimize potential contributions from actual electrical stimulation. Then, for each cluster, we ran a robust correlation (Rousselet & Pernet, 2012; Wilcox, 2012) across participants. For each participant, we considered the average regression coefficient corresponding to the proximity by direction interaction (from the original model so as to estimate it with minimal contamination from shocks) and regression coefficient corresponding to physical shock events.

### Plotting Parametric Effects as a Function of Proximity

Equation 1 allowed us to estimate the contributions of the seven main regressors to fMRI responses. Because of the parametric nature of the design, to illustrate responses in a more intuitive manner, we estimated responses separately for approach and retreat for a range of proximity values (Figure 8). To do so, the value of z scored proximity was varied (in the range of [−2, 1.5] and at the mean speed value), and the estimated regression coefficients were used to estimate the response at each value of proximity.

To provide an indication of variability of the fit across participants, we adopted the following approach. In the case of the proximity by direction interaction (Figures 8 and 11A), at each level of proximity, we calculated the difference between the estimated response for the approach and retreat conditions. We then calculated the standard error of the approach-minus-retreat difference across participants (at each value of proximity). We display the 95% confidence bands at each proximity value (note that because the intervals were based on differences between approach and retreat conditions, the same band widths are used for approach and retreat). An analogous procedure was used for the proximity by direction interaction of SCRs (Figure 4). The BST exhibited a proximity effect but no interaction. Therefore, in Figure 9 we computed error bands separately for approach and retreat based on the variability of estimated responses across participants as a function of proximity.

Figure 4.

Skin conductance response proximity by direction interaction. Estimated responses for a range of proximity values. To display estimated responses, we varied proximity and estimated the response based on the linear model for SCR (analogous to the model of Equation 1). The approach versus retreat difference was greater when circles were near compared with far. The confidence bands were obtained by considering within-subject differences (approach minus retreat); see Methods. A.U. = arbitrary units.

Figure 4.

Skin conductance response proximity by direction interaction. Estimated responses for a range of proximity values. To display estimated responses, we varied proximity and estimated the response based on the linear model for SCR (analogous to the model of Equation 1). The approach versus retreat difference was greater when circles were near compared with far. The confidence bands were obtained by considering within-subject differences (approach minus retreat); see Methods. A.U. = arbitrary units.

### Statistical Approach and p Values

The null hypothesis significance testing framework has come under increased scrutiny in recent years. In particular, the hard threshold of .05 has come under attack, with reasonable researchers calling for both stricter thresholds (Benjamin et al., 2018) or, conversely, for p values to be abandoned (McShane, Gal, Gelman, Robert, & Tackett, 2017). However, like McShane and colleagues, we do not consider a binary threshold to be satisfactory and believe that p values should be treated continuously. Accordingly, in select cases, we show p values and discuss findings that do not survive correction for multiple comparisons; in the context of Table 9, we discuss the general results of the BST given its important role in threat-related processing.

## RESULTS

Our paradigm allowed us to investigate the role played by threat proximity, direction, and speed and their interactions on SCRs and fMRI responses. Intuitively, interactions evaluated the extent to which factor combinations were relevant in explaining the data. For instance, the contrast of approach versus retreat (direction) was anticipated to depend on proximity (Figure 1C). Moreover, as proximity and speed varied continuously, their roles and their interactions were assessed parametrically.

Our design did not include a standard control condition (e.g., circles colliding but no shock administered), as often is the case in fMRI studies. Note, however, that our main goal was not to investigate the shock event itself but potential threat. Thus, approach and retreat can be viewed as paired conditions insofar as processes related to tracking the movement of the circles are concerned, for example. Furthermore, as stated in the preceding paragraph, an important focus of the research was to assess whether or not brain regions were sensitive to variable interactions, an approach that further helped reduce the contributions of non-threat related processing (see also Discussion).

### Skin Conductance Responses

Analysis of SCR data revealed that all three main variables had robust effects on responses (Table 1). In addition, we detected an interaction of proximity by direction; in this case, responses to approach versus retreat were sensitive to threat distance, such that the effect was larger when near versus far. To visualize this result, Figure 4 shows estimated SCRs for approach and retreat for a range of proximity values (because the circles moved continuously on the screen, Figure 4 used an approach similar to that of Figure 8 for plotting; see Methods). Finally, a three-way interaction between proximity, direction, and speed also survived correction for multiple comparisons.

Table 1.
SCR Results
Regressort(81)p
Proximity 4.57 .0000
Direction 9.37 .0000
Speed −4.20 .0001
Direction × Speed −0.92 .3602
Proximity × Direction 10.99 .0000
Proximity × Speed −2.43 .0175
Proximity × Direction × Speed −2.78 .0067
Regressort(81)p
Proximity 4.57 .0000
Direction 9.37 .0000
Speed −4.20 .0001
Direction × Speed −0.92 .3602
Proximity × Direction 10.99 .0000
Proximity × Speed −2.43 .0175
Proximity × Direction × Speed −2.78 .0067

Bonferroni correction for multiple comparisons: 0.05/7 = 0.0071.

### fMRI Voxelwise Analysis

Figures 5 and 6 (Tables 2 and 3) show the effects of proximity and direction (Table 4 shows the effect of speed). The main focus of this study was to investigate interactions between threat-related factors. Figure 7 (Table 5) shows interactions between proximity and direction; positive voxels (red) show effects when the contrast of approach versus retreat was greater during closer versus farther circles, and blue voxels indicate the opposite. Figure 8 shows estimated responses for approach and retreat for a range of proximity values, which aids in visualizing the parametric effects of proximity on the signals in the two regions (see Methods). For the right anterior insula (Figure 8A), when the circles were closer to each other, a larger approach versus retreat differential response was observed compared with when the circles were farther from each other. Responses for the right amygdala (Figure 8B) exhibited the opposite pattern as responses were larger for retreat compared with approach, and the contrast was enhanced when circles were closer compared with farther. Tables 6 and 7 show two-way interactions between direction and speed and between proximity and speed. Table 8 shows the three-way interaction of proximity, direction, and speed.

Figure 5.

Brain responses as a function of threat proximity. Clusters in red show regions with stronger responses for closer versus farther; clusters in blue show the reverse. Clusters were thresholded at a whole-brain corrected alpha of .05. SMA = supplementary motor area; FEF = frontal eye field; IPS = intraparietal sulcus; SFG = superior frontal gyrus; PCC = posterior cingulate cortex; vmPFC = ventromedial prefrontal cortex.

Figure 5.

Brain responses as a function of threat proximity. Clusters in red show regions with stronger responses for closer versus farther; clusters in blue show the reverse. Clusters were thresholded at a whole-brain corrected alpha of .05. SMA = supplementary motor area; FEF = frontal eye field; IPS = intraparietal sulcus; SFG = superior frontal gyrus; PCC = posterior cingulate cortex; vmPFC = ventromedial prefrontal cortex.

Figure 6.

Brain responses as a function of direction (approach vs. retreat). Clusters in red show regions with stronger responses for approach versus retreat; clusters in blue show the reverse. Clusters were thresholded at a whole-brain corrected alpha of .05. PAG = periaqueductal gray; SMA = supplementary motor area; FEF = frontal eye field; IPS = intraparietal sulcus; PreCG = precentral gyrus.

Figure 6.

Brain responses as a function of direction (approach vs. retreat). Clusters in red show regions with stronger responses for approach versus retreat; clusters in blue show the reverse. Clusters were thresholded at a whole-brain corrected alpha of .05. PAG = periaqueductal gray; SMA = supplementary motor area; FEF = frontal eye field; IPS = intraparietal sulcus; PreCG = precentral gyrus.

Table 2.
Clusters that Exhibited the Effect of Proximity in Voxelwise Analysis at Whole-brain Cluster-level Corrected Alpha of .05
kxyztCluster
12470 −14 −88 28 −13.47 Occipital cortex/cuneus/posterior cingulate cortex
1690 36 22 7.93 Right anterior/mid-insula
1489 −34 −92 −6 8.84 Left inferior/middle occipital gyrus
1453 28 −90 8.93 Right inferior/middle occipital gyrus
1188 64 −38 30 7.83 Right supramarginal/postcentral gyrus
1088 14 10 64 6.40 Right superior frontal gyrus
995 −16 −76 −34 7.84 Left cerebellum
869 −60 −46 42 6.58 Left supramarginal gyrus
796 −32 22 7.28 Left anterior insula
576 −2 46 −10 −5.91 Ventromedial prefrontal cortex
526 −18 7.49 Right/left thalamus
336 −22 26 46 −5.23 Left superior frontal gyrus
333 −12 −72 −44 5.89 Left cerebellum
209 44 −58 −30 5.88 Right cerebellum
138 22 32 48 −4.55 Right superior frontal gyrus
125 −4 −20 30 5.27 Right/left posterior cingulate cortex
118 −34 48 28 4.96 Left middle frontal gyrus
117 −62 −6 −14 −5.84 Left middle temporal gyrus
117 18 18 5.41 Right dorsolateral caudate
88 26 42 22 4.47 Right middle frontal gyrus
83 −26 −10 5.18 Left putamen
79 −12 −54 66 −5.71 Left precuneus/superior parietal lobule
72 32 48 4.88 Medial superior frontal gyrus
70 58 −30 −6 4.36 Right superior/middle temporal gyrus
68 56 −4 −18 −4.98 Right middle temporal gyrus
66 −26 −24 54 −4.83 Left precentral gyrus
65 −34 −6 50 4.57 Left middle frontal gyrus
61 42 −12 48 −4.85 Right precentral gyrus
55 −42 −24 18 −4.11 Left posterior insula
42 32 40 −10 −5.02 Right lateral orbitofrontal cortex
40 −12 −2 66 4.73 Left superior frontal gyrus
36 −58 −22 48 −4.09 Left postcentral gyrus
34 36 −70 −10 5.04 Right inferior temporal gyrus
34 12 −10 −5.17 Subcollosal area
34 −4 −68 50 −4.38 Precuneus
31 20 −14 −24 −5.47 Right hippocampus/amygdala
25 −22 −8 −22 −4.36 Left hippocampus/amygdala
24 −4 −24 54 −4.52 Left paracentral lobule
20 38 −12 −8 4.36 Right planum polare
20 54 −26 10 −4.20 Right tranverse temporal gyrus
17 12 12 −2 5.02 Right ventral caudate
17 −52 −28 10 −3.91 Left tranverse temporal gyrus
17 −16 64 14 −4.22 Right superior frontopolar gyrus
17 66 −4 18 −4.78 Right precentral gyrus
17 14 −54 68 −4.57 Right precuneus/superior parietal lobule
16 −26 62 −12 4.42 Left frontomarginal gyrus
15 12 −74 −40 5.49 Right cerebellum
14 −56 −50 −10 −4.11 Left superior/middle temporal gyrus
13 −18 −34 72 −4.34 Left postcentral gyrus
kxyztCluster
12470 −14 −88 28 −13.47 Occipital cortex/cuneus/posterior cingulate cortex
1690 36 22 7.93 Right anterior/mid-insula
1489 −34 −92 −6 8.84 Left inferior/middle occipital gyrus
1453 28 −90 8.93 Right inferior/middle occipital gyrus
1188 64 −38 30 7.83 Right supramarginal/postcentral gyrus
1088 14 10 64 6.40 Right superior frontal gyrus
995 −16 −76 −34 7.84 Left cerebellum
869 −60 −46 42 6.58 Left supramarginal gyrus
796 −32 22 7.28 Left anterior insula
576 −2 46 −10 −5.91 Ventromedial prefrontal cortex
526 −18 7.49 Right/left thalamus
336 −22 26 46 −5.23 Left superior frontal gyrus
333 −12 −72 −44 5.89 Left cerebellum
209 44 −58 −30 5.88 Right cerebellum
138 22 32 48 −4.55 Right superior frontal gyrus
125 −4 −20 30 5.27 Right/left posterior cingulate cortex
118 −34 48 28 4.96 Left middle frontal gyrus
117 −62 −6 −14 −5.84 Left middle temporal gyrus
117 18 18 5.41 Right dorsolateral caudate
88 26 42 22 4.47 Right middle frontal gyrus
83 −26 −10 5.18 Left putamen
79 −12 −54 66 −5.71 Left precuneus/superior parietal lobule
72 32 48 4.88 Medial superior frontal gyrus
70 58 −30 −6 4.36 Right superior/middle temporal gyrus
68 56 −4 −18 −4.98 Right middle temporal gyrus
66 −26 −24 54 −4.83 Left precentral gyrus
65 −34 −6 50 4.57 Left middle frontal gyrus
61 42 −12 48 −4.85 Right precentral gyrus
55 −42 −24 18 −4.11 Left posterior insula
42 32 40 −10 −5.02 Right lateral orbitofrontal cortex
40 −12 −2 66 4.73 Left superior frontal gyrus
36 −58 −22 48 −4.09 Left postcentral gyrus
34 36 −70 −10 5.04 Right inferior temporal gyrus
34 12 −10 −5.17 Subcollosal area
34 −4 −68 50 −4.38 Precuneus
31 20 −14 −24 −5.47 Right hippocampus/amygdala
25 −22 −8 −22 −4.36 Left hippocampus/amygdala
24 −4 −24 54 −4.52 Left paracentral lobule
20 38 −12 −8 4.36 Right planum polare
20 54 −26 10 −4.20 Right tranverse temporal gyrus
17 12 12 −2 5.02 Right ventral caudate
17 −52 −28 10 −3.91 Left tranverse temporal gyrus
17 −16 64 14 −4.22 Right superior frontopolar gyrus
17 66 −4 18 −4.78 Right precentral gyrus
17 14 −54 68 −4.57 Right precuneus/superior parietal lobule
16 −26 62 −12 4.42 Left frontomarginal gyrus
15 12 −74 −40 5.49 Right cerebellum
14 −56 −50 −10 −4.11 Left superior/middle temporal gyrus
13 −18 −34 72 −4.34 Left postcentral gyrus

Peak MNI coordinates, t(84) values, and cluster size (k) refer to number of 2.0 × 2.0 × 2.0 mm3 voxels. Peak coordinates are presented for completeness and potential meta-analysis; with cluster-based thresholding, it is not possible to conclude that all the reported peaks were activated (see Woo, Krishnan, & Wager, 2014).

Table 3.
Clusters that Exhibited the Effect of Direction in Voxelwise Analysis at Whole-brain Cluster-level Corrected Alpha of .05
kxyztCluster
761 −36 −44 56 7.77 Left inferior parietal cortex
633 34 −50 60 6.67 Right inferior parietal cortex
572 −46 −62 12 7.29 Left superior temporal gyrus
565 −24 −8 52 8.40 Left frontal eye field
546 48 −60 10 6.78 Right superior temporal gyrus
472 36 −4 50 6.62 Right frontal eye field
277 26 −98 −4 −6.22 Right superior occipital gyrus
269 −28 −72 26 7.02 Left parieto-occipitalis
246 26 −68 −4 −6.49 Right lingual gyrus
230 36 28 7.45 Right anterior insula
205 14 −86 24 −6.17 Right parieto-occipitalis (posterior)
151 −26 −98 −4 −6.63 Left superior occipital gyrus
144 −56 38 6.91 Left precentral gyrus
134 −34 18 6.73 Left anterior insula
130 54 34 5.89 Right precentral gyrus
108 56 −44 20 4.73 Right parietal operculum
106 −16 −86 22 −5.48 Left parieto-occipitalis (posterior)
105 −12 −94 16 −5.63 Left parieto-occipitalis (posterior)
96 32 −72 32 5.32 Right parieto-occipitalis
90 −30 −28 52 −4.79 Left postcentral gyrus
82 −12 −76 −8 −5.77 Left lingual gyrus
63 −28 −58 −8 −4.87 Left lingual/fusiform gyrus
60 14 −80 −4.54 Right occipital gyrus
43 12 −70 −20 5.39 Right cerebellum
38 −28 −10 4.71 Right periaqueductal gray
37 −8 −74 −42 5.30 Left cerebellum
37 −42 −74 −8 4.35 Left inferior temporal gyrus
35 46 20 26 4.88 Right inferior/middle frontal gyrus
26 −40 −54 −18 4.46 Left fusiform gyrus
26 −12 −46 52 5.04 Left paracentral lobule
24 −12 −22 40 4.68 Left posterior cingulate cortex
23 44 −46 −14 4.93 Right fusiform gyrus
22 10 −74 −38 6.09 Right cerebellum
22 −54 −32 4.59 Cerebellum
22 −20 20 48 −4.80 Left superior frontal gyrus
21 20 −24 66 −4.74 Right postcentral gyrus
19 26 −10 4.60 Right putamen
18 −18 −72 −22 4.60 Left cerebellum
18 12 70 4.81 Right superior frontal gyrus
16 36 −14 44 −4.56 Right precentral gyrus
13 42 42 −2 −4.33 Right inferior frontal gyrus
13 −86 18 −4.80 Right parieto-occipitalis (posterior)
kxyztCluster
761 −36 −44 56 7.77 Left inferior parietal cortex
633 34 −50 60 6.67 Right inferior parietal cortex
572 −46 −62 12 7.29 Left superior temporal gyrus
565 −24 −8 52 8.40 Left frontal eye field
546 48 −60 10 6.78 Right superior temporal gyrus
472 36 −4 50 6.62 Right frontal eye field
277 26 −98 −4 −6.22 Right superior occipital gyrus
269 −28 −72 26 7.02 Left parieto-occipitalis
246 26 −68 −4 −6.49 Right lingual gyrus
230 36 28 7.45 Right anterior insula
205 14 −86 24 −6.17 Right parieto-occipitalis (posterior)
151 −26 −98 −4 −6.63 Left superior occipital gyrus
144 −56 38 6.91 Left precentral gyrus
134 −34 18 6.73 Left anterior insula
130 54 34 5.89 Right precentral gyrus
108 56 −44 20 4.73 Right parietal operculum
106 −16 −86 22 −5.48 Left parieto-occipitalis (posterior)
105 −12 −94 16 −5.63 Left parieto-occipitalis (posterior)
96 32 −72 32 5.32 Right parieto-occipitalis
90 −30 −28 52 −4.79 Left postcentral gyrus
82 −12 −76 −8 −5.77 Left lingual gyrus
63 −28 −58 −8 −4.87 Left lingual/fusiform gyrus
60 14 −80 −4.54 Right occipital gyrus
43 12 −70 −20 5.39 Right cerebellum
38 −28 −10 4.71 Right periaqueductal gray
37 −8 −74 −42 5.30 Left cerebellum
37 −42 −74 −8 4.35 Left inferior temporal gyrus
35 46 20 26 4.88 Right inferior/middle frontal gyrus
26 −40 −54 −18 4.46 Left fusiform gyrus
26 −12 −46 52 5.04 Left paracentral lobule
24 −12 −22 40 4.68 Left posterior cingulate cortex
23 44 −46 −14 4.93 Right fusiform gyrus
22 10 −74 −38 6.09 Right cerebellum
22 −54 −32 4.59 Cerebellum
22 −20 20 48 −4.80 Left superior frontal gyrus
21 20 −24 66 −4.74 Right postcentral gyrus
19 26 −10 4.60 Right putamen
18 −18 −72 −22 4.60 Left cerebellum
18 12 70 4.81 Right superior frontal gyrus
16 36 −14 44 −4.56 Right precentral gyrus
13 42 42 −2 −4.33 Right inferior frontal gyrus
13 −86 18 −4.80 Right parieto-occipitalis (posterior)

Peak MNI coordinates, t(84) values, and cluster size (k) refer to number of 2.0 × 2.0 × 2.0 mm3 voxels. Peak coordinates are presented for completeness and potential meta-analysis; with cluster-based thresholding, it is not possible to conclude that all the reported peaks were activated (see Woo et al., 2014).

Table 4.
Clusters that Exhibited the Effect of Speed in Voxelwise Analysis at Whole-brain Cluster-level Corrected Alpha of .05
kxyztCluster
4402 −46 −74 9.74 Left inferior/middle temporal gyrus/fusiform gyrus
3332 46 −68 9.38 Right inferior/middle temporal gyrus/fusiform gyrus
339 −26 −56 52 5.70 Left intraparietal sulcus
283 −32 −4 46 5.68 Left frontal eye field
235 26 −50 52 6.23 Right intraparietal sulcus
178 −32 24 6.92 Left anterior insula
120 36 −4 50 5.35 Right frontal eye field
114 36 24 5.42 Right anterior insula
112 −6 50 5.52 Left mid-cingulate cortex/supplementary motor area
110 −12 −24 40 6.08 Left posterior cingulate cortex
102 −50 36 6.63 Left precentral gyrus
63 10 20 36 5.22 Right mid-cingulate cortex
50 54 −44 20 4.98 Right parietal operculum
45 56 4.50 Right supplementary motor area
42 50 34 5.36 Right precentral gyrus
35 12 −94 20 4.48 Right parieto-occipitalis (posterior)
33 18 64 4.66 Right superior frontal gyrus
29 20 −74 40 4.81 Right parieto-occipitalis
28 −34 −46 −20 5.21 Left fusiform gyrus
25 −12 −74 12 4.56 Left precuneus/occipital gyrus
24 −14 −30 −2 5.14 Left ventral thalamus
19 −14 −46 48 4.43 Left superior parietal lobule
17 −36 −12 −6 4.36 Left postcentral insular cortex
16 24 −70 10 4.07 Right precuneus/occipital gyrus
15 12 −20 40 4.00 Right posterior cingulate cortex
14 −8 −72 −38 5.52 Left cerebellum
13 −8 26 30 4.10 Left mid-cingulate cortex
13 36 34 4.26 Right precentral gyrus
13 −50 −26 36 4.50 Left supramarginal gyrus
kxyztCluster
4402 −46 −74 9.74 Left inferior/middle temporal gyrus/fusiform gyrus
3332 46 −68 9.38 Right inferior/middle temporal gyrus/fusiform gyrus
339 −26 −56 52 5.70 Left intraparietal sulcus
283 −32 −4 46 5.68 Left frontal eye field
235 26 −50 52 6.23 Right intraparietal sulcus
178 −32 24 6.92 Left anterior insula
120 36 −4 50 5.35 Right frontal eye field
114 36 24 5.42 Right anterior insula
112 −6 50 5.52 Left mid-cingulate cortex/supplementary motor area
110 −12 −24 40 6.08 Left posterior cingulate cortex
102 −50 36 6.63 Left precentral gyrus
63 10 20 36 5.22 Right mid-cingulate cortex
50 54 −44 20 4.98 Right parietal operculum
45 56 4.50 Right supplementary motor area
42 50 34 5.36 Right precentral gyrus
35 12 −94 20 4.48 Right parieto-occipitalis (posterior)
33 18 64 4.66 Right superior frontal gyrus
29 20 −74 40 4.81 Right parieto-occipitalis
28 −34 −46 −20 5.21 Left fusiform gyrus
25 −12 −74 12 4.56 Left precuneus/occipital gyrus
24 −14 −30 −2 5.14 Left ventral thalamus
19 −14 −46 48 4.43 Left superior parietal lobule
17 −36 −12 −6 4.36 Left postcentral insular cortex
16 24 −70 10 4.07 Right precuneus/occipital gyrus
15 12 −20 40 4.00 Right posterior cingulate cortex
14 −8 −72 −38 5.52 Left cerebellum
13 −8 26 30 4.10 Left mid-cingulate cortex
13 36 34 4.26 Right precentral gyrus
13 −50 −26 36 4.50 Left supramarginal gyrus

Peak MNI coordinates, t(84) values, and cluster size (k) refer to number of 2.0 × 2.0 × 2.0 mm3 voxels. Peak coordinates are presented for completeness and potential meta-analysis; with cluster-based thresholding, it is not possible to conclude that all the reported peaks were activated (see Woo et al., 2014).

Figure 7.

Brain responses exhibiting a proximity by direction (approach vs. retreat) interaction in areas of interest. Clusters in red show regions with approach versus retreat responses greater when closer versus farther; clusters in blue show the reserve pattern. Clusters were thresholded at a whole-brain corrected alpha of .05. FEF = frontal eye field; PreCG = precentral gyrus.

Figure 7.

Brain responses exhibiting a proximity by direction (approach vs. retreat) interaction in areas of interest. Clusters in red show regions with approach versus retreat responses greater when closer versus farther; clusters in blue show the reserve pattern. Clusters were thresholded at a whole-brain corrected alpha of .05. FEF = frontal eye field; PreCG = precentral gyrus.

Table 5.
Clusters that Exhibited the Proximity × Direction Interaction in Voxelwise Analysis at Whole-brain Cluster-level Corrected Alpha of .05
kxyztCluster
528 10 −74 −2 −9.42 Right occipital cortex
398 −10 −96 12 −7.09 Left occipital gyrus
358 28 −96 −5.98 Right occipital gyrus
263 −20 −12 60 6.26 Left frontal eye field
243 16 −86 26 −7.16 Right occipital gyrus
215 26 −2 58 5.55 Right frontal eye field
159 −54 −32 −2 −5.16 Left superior temporal gyrus
148 −16 −86 22 −5.94 Left occipital gyrus
138 −28 −98 −5.24 Left inferior occipital gyrus
109 −54 38 6.72 Left precentral gyrus
109 50 −32 58 −4.88 Right postcentral gyrus
105 16 −92 18 −5.51 Right occipital gyrus
99 66 −16 20 −4.95 Right supramarginal gyrus
98 −30 −28 52 −4.83 Left precentral gyrus
74 34 28 5.63 Right anterior insula
73 40 −60 50 −4.61 Right angular gyrus
54 22 −32 72 −5.81 Right postcentral gyrus
52 20 −60 12 −4.85 Right occipital gyrus
35 −44 62 −4.42 Right superior postcentral sulcus
34 18 −74 24 −5.00 Right posterior angular gyrus
31 −18 −66 −5.02 Left occipital gyrus
31 −36 −16 16 −4.80 Left posterior insula
29 54 36 5.11 Right precentral gyrus
27 −54 −22 54 −4.95 Left postcentral gyrus
25 −26 −70 −28 −4.45 Left cerebellum
24 −60 −32 16 −4.51 Left supramarginal gyrus
24 38 −16 40 −5.01 Right postcentral gyrus
23 56 −58 −6 −4.01 Right inferior temporal gyrus
22 28 −10 −20 −4.16 Right amygdala
22 −50 22 −8 −4.66 Left inferior temporal gyrus
21 −36 −18 44 −5.06 Left postcentral gyrus
20 −54 24 18 −4.34 Left inferior frontal gyrus
20 −50 −20 18 −4.87 Left parietal operculum
20 −2 22 54 −4.24 Left paracentral lobule
19 32 −76 −38 −4.32 Right cerebellum
18 58 −58 28 −4.07 Right supramarginal gyrus
17 44 −70 −20 −3.88 Right inferior temporal gyrus
15 54 20 −4 −4.59 Right inferior temporal gyrus
15 −8 −76 16 −4.19 Left precuneus
15 32 −36 50 4.21 Right postcentral gyrus
14 −52 22 24 −3.98 Left middle frontal gyrus
13 −36 −60 −42 −3.70 Left cerebellum
kxyztCluster
528 10 −74 −2 −9.42 Right occipital cortex
398 −10 −96 12 −7.09 Left occipital gyrus
358 28 −96 −5.98 Right occipital gyrus
263 −20 −12 60 6.26 Left frontal eye field
243 16 −86 26 −7.16 Right occipital gyrus
215 26 −2 58 5.55 Right frontal eye field
159 −54 −32 −2 −5.16 Left superior temporal gyrus
148 −16 −86 22 −5.94 Left occipital gyrus
138 −28 −98 −5.24 Left inferior occipital gyrus
109 −54 38 6.72 Left precentral gyrus
109 50 −32 58 −4.88 Right postcentral gyrus
105 16 −92 18 −5.51 Right occipital gyrus
99 66 −16 20 −4.95 Right supramarginal gyrus
98 −30 −28 52 −4.83 Left precentral gyrus
74 34 28 5.63 Right anterior insula
73 40 −60 50 −4.61 Right angular gyrus
54 22 −32 72 −5.81 Right postcentral gyrus
52 20 −60 12 −4.85 Right occipital gyrus
35 −44 62 −4.42 Right superior postcentral sulcus
34 18 −74 24 −5.00 Right posterior angular gyrus
31 −18 −66 −5.02 Left occipital gyrus
31 −36 −16 16 −4.80 Left posterior insula
29 54 36 5.11 Right precentral gyrus
27 −54 −22 54 −4.95 Left postcentral gyrus
25 −26 −70 −28 −4.45 Left cerebellum
24 −60 −32 16 −4.51 Left supramarginal gyrus
24 38 −16 40 −5.01 Right postcentral gyrus
23 56 −58 −6 −4.01 Right inferior temporal gyrus
22 28 −10 −20 −4.16 Right amygdala
22 −50 22 −8 −4.66 Left inferior temporal gyrus
21 −36 −18 44 −5.06 Left postcentral gyrus
20 −54 24 18 −4.34 Left inferior frontal gyrus
20 −50 −20 18 −4.87 Left parietal operculum
20 −2 22 54 −4.24 Left paracentral lobule
19 32 −76 −38 −4.32 Right cerebellum
18 58 −58 28 −4.07 Right supramarginal gyrus
17 44 −70 −20 −3.88 Right inferior temporal gyrus
15 54 20 −4 −4.59 Right inferior temporal gyrus
15 −8 −76 16 −4.19 Left precuneus
15 32 −36 50 4.21 Right postcentral gyrus
14 −52 22 24 −3.98 Left middle frontal gyrus
13 −36 −60 −42 −3.70 Left cerebellum

Peak MNI coordinates, t(84) values, and cluster size (k) refer to number of 2.0 × 2.0 × 2.0 mm3 voxels. Peak coordinates are presented for completeness and potential meta-analysis; with cluster-based thresholding, it is not possible to conclude that all the reported peaks were activated (see Woo et al., 2014).

Figure 8.

Proximity by direction (approach vs. retreat) interaction. Estimated responses for a range of proximity values. (A) For the right anterior insula, activity increased as a function of proximity for both approach and retreat, but more steeply for the former. (B) For the right amygdala, activity decreased as a function of proximity during approach, but changed little during retreat. The confidence bands were obtained by considering within-subject differences (approach minus retreat); see Methods. A.U. = arbitrary units.

Figure 8.

Proximity by direction (approach vs. retreat) interaction. Estimated responses for a range of proximity values. (A) For the right anterior insula, activity increased as a function of proximity for both approach and retreat, but more steeply for the former. (B) For the right amygdala, activity decreased as a function of proximity during approach, but changed little during retreat. The confidence bands were obtained by considering within-subject differences (approach minus retreat); see Methods. A.U. = arbitrary units.

Table 6.
Clusters that Exhibited the Direction × Speed Interaction in Voxelwise Analysis at Whole-brain Cluster-level Corrected Alpha of .05
kxyztCluster
125 −46 −74 4.80 Left middle temporal gyrus
105 −30 −100 −2 −5.40 Left superior occipital gyrus
95 34 −94 −4.82 Right middle occipital gyrus
95 −32 −26 58 −4.34 Left precentral gyrus
55 34 26 4.44 Right anterior insula
46 16 −26 68 −4.54 Right precentral gyrus
30 −8 −22 62 −4.53 Left paracentral lobule
28 −30 −46 42 4.33 Left inferior parietal cortex
22 −24 −78 30 5.11 Left parieto-occipitalis
22 −8 −72 34 4.06 Left precuenus
22 −40 −32 42 4.53 Left inferior postcentral sulcus
21 −76 −2 −4.71 Right lingual gyrus
18 −30 26 −4 4.68 Left anterior insula
18 16 34 56 −4.58 Right superior frontal gyrus
17 −26 −54 −8 −3.94 Left parahippocampal gyrus
17 66 −2 16 −4.99 Right precentral gyrus
15 32 4.34 Mid-cingulate cortex
15 40 −2 48 4.21 Right precentral gyrus
15 30 −26 60 −4.11 Right precentral gyrus
14 −30 26 4.52 Left anterior insula
14 −4 −24 56 −4.55 Left paracentral lobule
13 46 48 −8 −4.55 Right inferior frontal gyrus (orbital)
13 20 46 36 −4.54 Right superior frontal gyrus
13 −4 −78 42 3.94 Left precuneus
kxyztCluster
125 −46 −74 4.80 Left middle temporal gyrus
105 −30 −100 −2 −5.40 Left superior occipital gyrus
95 34 −94 −4.82 Right middle occipital gyrus
95 −32 −26 58 −4.34 Left precentral gyrus
55 34 26 4.44 Right anterior insula
46 16 −26 68 −4.54 Right precentral gyrus
30 −8 −22 62 −4.53 Left paracentral lobule
28 −30 −46 42 4.33 Left inferior parietal cortex
22 −24 −78 30 5.11 Left parieto-occipitalis
22 −8 −72 34 4.06 Left precuenus
22 −40 −32 42 4.53 Left inferior postcentral sulcus
21 −76 −2 −4.71 Right lingual gyrus
18 −30 26 −4 4.68 Left anterior insula
18 16 34 56 −4.58 Right superior frontal gyrus
17 −26 −54 −8 −3.94 Left parahippocampal gyrus
17 66 −2 16 −4.99 Right precentral gyrus
15 32 4.34 Mid-cingulate cortex
15 40 −2 48 4.21 Right precentral gyrus
15 30 −26 60 −4.11 Right precentral gyrus
14 −30 26 4.52 Left anterior insula
14 −4 −24 56 −4.55 Left paracentral lobule
13 46 48 −8 −4.55 Right inferior frontal gyrus (orbital)
13 20 46 36 −4.54 Right superior frontal gyrus
13 −4 −78 42 3.94 Left precuneus

Peak MNI coordinates, t(84) values, and cluster size (k) refer to number of 2.0 × 2.0 × 2.0 mm3 voxels. Peak coordinates are presented for completeness and potential meta-analysis; with cluster-based thresholding, it is not possible to conclude that all the reported peaks were activated (see Woo et al., 2014).

Table 7.
Clusters that Exhibited the Proximity × Speed Interaction in Voxelwise Analysis at Whole-brain Cluster-level Corrected Alpha of .05
kxyztCluster
116 −2 56 5.58 Supplementary motor area
50 −6 16 34 4.68 Left mid-cingulate cortex
46 48 −66 −5.23 Right middle temporal gyrus
39 26 −90 −5.35 Right occipital gyrus
37 −30 28 5.25 Left anterior insula
33 30 −74 38 −4.46 Right parieto-occipitalis
23 40 −80 20 −4.38 Right occipital gyrus
22 22 38 4.65 Right mid-cingulate cortex
21 −62 4.19 Left precentral gyrus
19 18 −66 52 −4.47 Right superior parietal lobule
15 46 −58 −6 −4.49 Right inferior temporal gyrus
13 −24 −96 −4.09 Left middle occipital gyrus
13 12 −72 4.15 Right occipital gyrus
13 14 −54 64 −3.85 Right superior parietal lobule
kxyztCluster
116 −2 56 5.58 Supplementary motor area
50 −6 16 34 4.68 Left mid-cingulate cortex
46 48 −66 −5.23 Right middle temporal gyrus
39 26 −90 −5.35 Right occipital gyrus
37 −30 28 5.25 Left anterior insula
33 30 −74 38 −4.46 Right parieto-occipitalis
23 40 −80 20 −4.38 Right occipital gyrus
22 22 38 4.65 Right mid-cingulate cortex
21 −62 4.19 Left precentral gyrus
19 18 −66 52 −4.47 Right superior parietal lobule
15 46 −58 −6 −4.49 Right inferior temporal gyrus
13 −24 −96 −4.09 Left middle occipital gyrus
13 12 −72 4.15 Right occipital gyrus
13 14 −54 64 −3.85 Right superior parietal lobule

Peak MNI coordinates, t(84) values, and cluster size (k) refer to number of 2.0 × 2.0 × 2.0 mm3 voxels. Peak coordinates are presented for completeness and potential meta-analysis; with cluster-based thresholding, it is not possible to conclude that all the reported peaks were activated (see Woo et al., 2014).

Table 8.
Clusters that Exhibited the Proximity × Direction × Speed Interaction in Voxelwise Analysis at Whole-brain Cluster-level Corrected Alpha of .05
kxyztCluster
17 −54 −24 38 3.98 Left central sulcus
17 24 30 52 −5.02 Right superior frontal gyrus
kxyztCluster
17 −54 −24 38 3.98 Left central sulcus
17 24 30 52 −5.02 Right superior frontal gyrus

Peak MNI coordinates, t(84) values, and cluster size (k) refer to number of 2.0 × 2.0 × 2.0 mm3 voxels. Peak coordinates are presented for completeness and potential meta-analysis; with cluster-based thresholding, it is not possible to conclude that all the reported peaks were activated (see Woo et al., 2014).

### BST ROI Analysis

Given that the BST is a rather small region that is involved in threat-related processing, we ran a focused ROI analysis using anatomically defined left/right BST masks and unsmoothed data to minimize the influence of signals from surrounding structures. We observed a robust effect of threat proximity in the right BST (and weak evidence in the left BST), with stronger responses when circles were closer than farther (Figure 9A; Table 9). For the right BST, some evidence for proximity by speed interaction was seen.

Figure 9.

Proximity effect in the BST ROI analysis. Estimated responses for a range of proximity values. Activity increased as a function of proximity for both approach and retreat. The confidence bands were obtained by considering variability during approach and retreat, separately; see Methods. A.U. = arbitrary units.

Figure 9.

Proximity effect in the BST ROI analysis. Estimated responses for a range of proximity values. Activity increased as a function of proximity for both approach and retreat. The confidence bands were obtained by considering variability during approach and retreat, separately; see Methods. A.U. = arbitrary units.

Table 9.
BST ROI Analysis Results
RegressorLeft BSTRight BST
t(84)pt(84)p
Proximity 1.91 .0591 4.17 .0000
Direction 1.29 .1997 0.64 .0000
Speed 0.83 .4099 1.78 .0001
Direction × Speed −0.91 .3664 −0.26 .3602
Proximity × Direction −1.63 .1066 −0.35 .0000
Proximity × Speed −0.08 .9353 2.60 .0175
Proximity × Direction × Speed 0.00 .9986 −1.69 .0067
RegressorLeft BSTRight BST
t(84)pt(84)p
Proximity 1.91 .0591 4.17 .0000
Direction 1.29 .1997 0.64 .0000
Speed 0.83 .4099 1.78 .0001
Direction × Speed −0.91 .3664 −0.26 .3602
Proximity × Direction −1.63 .1066 −0.35 .0000
Proximity × Speed −0.08 .9353 2.60 .0175
Proximity × Direction × Speed 0.00 .9986 −1.69 .0067

Bonferroni correction for multiple comparisons: 0.05/7 = 0.0071.

### Relationship between SCR Responses and Brain Activity

We evaluated the linear relationship between SCR and fMRI by running a robust correlation analysis (across participants). Because multiple aspects of both the SCR and fMRI data could be probed (simple effects and interactions), we chose to focus the interrogation on the proximity by direction interaction. Thus, for both SCR and fMRI, the strength of the two-way interaction was considered for the analysis (as given by the regression coefficient in Equation 1). To minimize the problem of multiple statistical comparisons, for this analysis, we focused on clusters exhibiting a two-way interaction in the right anterior insula and the right amygdala, regions that feature in most models of threat processing. We did not detect a relationship between SCR and fMRI responses in either the right anterior insula, r(77) = .07, p = .550, or the right amygdala, r(75) = −.04, p = .697.

### Relationship between Anticipatory Activity and Physical Shock Responses

Our interpretation of the proximity by direction interaction was that it reflected, at least in part, threat-related processing, especially in brain regions important for this type of processing, such as the anterior insula. In an exploratory analysis, we tested if the strength of this interaction effect was associated (across participants) with the strength of responses evoked by physical shock. For the right anterior insula cluster that exhibited a proximity by direction interaction, we detected a positive linear relationship between the two measures, r(80) = .33, p = .002 (Figure 10). Given the importance of the amygdala in threat processing, we also tested the relationship in the right amygdala (also considering the cluster that exhibited a proximity by direction interaction), but no effect was detected, r(80) = −.02, p = .888.

Figure 10.

Relationship between anticipatory activity and physical shock responses in the right anterior insula. For the anticipatory activity, the proximity by direction interaction was considered for the analysis. Data points correspond to participants (red points indicate outliers deemed based on the robust correlation algorithm). A.U. = arbitrary units.

Figure 10.

Relationship between anticipatory activity and physical shock responses in the right anterior insula. For the anticipatory activity, the proximity by direction interaction was considered for the analysis. Data points correspond to participants (red points indicate outliers deemed based on the robust correlation algorithm). A.U. = arbitrary units.

### Individual Differences in State and Trait Anxiety

Linear relationships between state/trait anxiety and SCR or, separately, fMRI interactions of proximity and direction in the right anterior insula were not detected (all rs < .1 in absolute value). We detected a modest positive relationship between state anxiety and fMRI interactions of proximity and direction in the right amygdala (in the cluster that exhibited a proximity by direction interaction; state: r(77) = .2107, p = .0544; trait: r(79) = .0769, p = .4870). Given the multiple tests involved here, we do not believe these findings are noteworthy.

### Exploratory Analyses: PAG Responses

To visualize the responses of the PAG, we plotted estimated responses (Figure 11A), as done above for the right anterior insula, right amygdala, and right BST. To do so, we used the cluster (38 voxels) that exhibited the direction effect (approach vs. retreat) previously reported (Figure 6). Upon plotting, we discerned an effect of proximity for the approach condition, but not for retreat, consistent with a proximity by direction interaction (which was not detected in the voxelwise analysis). Given the importance of the PAG in the orchestration of defensive responses in the face of threat (Pessoa, 2016; Bandler & Shipley, 1994), we performed an additional exploratory analysis in this region. First, we generated a representative time series for the PAG by averaging the time series of the voxels within the cluster (based on the voxelwise effect of direction) and then evaluated the full model (Equation 1). As shown in Table 10, a robust proximity by direction interaction was detected (note that the interaction effects were nearly independent from the selection criterion, which was based on direction; the correlation between the interaction and direction was −.14). Given this result, we inspected again the results at the voxelwise level and observed some voxels that exhibited such an interaction, but too few to survive cluster thresholding.

Figure 11.

Exploratory analysis of the PAG. (A) Estimated responses for a range of proximity values. During approach, activity increased as a function of proximity; activity changed little during retreat periods. The confidence bands were obtained by considering within-subject differences (approach minus retreat); see Methods. A.U. = arbitrary units. (B) Contour plots show estimated responses for different combinations of proximity and speed during approach and retreat periods. Arrows point in the direction of signal increase. During approach, both proximity and speed simultaneously influenced responses, which increased when the circles were closer and speed was higher. A.U. = arbitrary units.

Figure 11.

Exploratory analysis of the PAG. (A) Estimated responses for a range of proximity values. During approach, activity increased as a function of proximity; activity changed little during retreat periods. The confidence bands were obtained by considering within-subject differences (approach minus retreat); see Methods. A.U. = arbitrary units. (B) Contour plots show estimated responses for different combinations of proximity and speed during approach and retreat periods. Arrows point in the direction of signal increase. During approach, both proximity and speed simultaneously influenced responses, which increased when the circles were closer and speed was higher. A.U. = arbitrary units.

Table 10.
Exploratory Right PAG ROI Analysis Results
Regressort(84)p
Proximity 2.33 .0220
Direction 6.75 .0000
Speed 2.69 .0087
Direction × Speed 1.92 .0588
Proximity × Direction 3.89 .0002
Proximity × Speed 1.17 .2455
Proximity × Direction × Speed 2.86 .0054
Regressort(84)p
Proximity 2.33 .0220
Direction 6.75 .0000
Speed 2.69 .0087
Direction × Speed 1.92 .0588
Proximity × Direction 3.89 .0002
Proximity × Speed 1.17 .2455
Proximity × Direction × Speed 2.86 .0054

Bonferroni correction for multiple comparisons: 0.05/7 = 0.0071.

Notably, we also observed a robust three-way interaction. As the three factors simultaneously affected PAG responses, the finding can be visualized via a contour plot (Figure 11B). During approach periods, when proximity increased (circles moved closer to each other), stronger responses were observed as speed increased from slower to faster (compare the top right vs. bottom left quadrants).

### Exploratory Analyses: Potential Nonlinear Effects of Proximity

The regression model we used (Equation 1) makes the assumption that the effect of proximity is linear. In additional exploratory analyses, we investigated potential nonlinear effects of proximity on brain activity. To do so, we inspected the pattern of the residuals as a function of proximity in the right anterior insula, right amygdala, right BST, and right PAG. For example, Figure 12 shows the residuals when using Equation 1 for the right anterior insula. Based on the pattern of residuals, the linear modeling approach adopted here appears to be reasonable in the context of our experiment.

Figure 12.

Exploratory analysis of potential nonlinear effects of proximity. The residuals from the model fit are plotted as a function of proximity. No appreciable lack of fit is evident. To plot residuals for all participants, they were first studentized (jitter as a function of proximity was also used to reduce overlap).

Figure 12.

Exploratory analysis of potential nonlinear effects of proximity. The residuals from the model fit are plotted as a function of proximity. No appreciable lack of fit is evident. To plot residuals for all participants, they were first studentized (jitter as a function of proximity was also used to reduce overlap).

## DISCUSSION

In this study, we investigated the role of threat-related factors and their temporally evolving interactions. Our findings support the view that threat processing is context sensitive and dynamic (Mobbs et al., 2015; Kavaliers & Choleris, 2001; Blanchard & Blanchard, 1990; Fanselow & Lester, 1988). In some brain regions, signal fluctuations were sensitive to continuous manipulations of proximity and speed indicating that threat processing is dynamic. Importantly, whereas some brain regions tracked individual threat-related factors (proximity, direction, or speed), others were also sensitive to combinations of these variables revealing the context-sensitive nature of threat processing. In this section, we will focus the discussion on a few of the brain regions that have been most heavily implicated in threat-related processing in the literature, specifically the anterior insula, amygdala, BST, and PAG.

To investigate how threat-related factors influence physiological arousal during dynamic threat anticipation, we recorded SCR during scanning. We observed robust effects of proximity and direction, with larger responses during near versus far and approach versus retreat, respectively. Of note, we observed a robust proximity by direction interaction, where responses to threat direction (approach vs. retreat) were enhanced when the circles were “near” compared with “far,” suggesting that the influence of dynamic threat anticipation on physiological arousal was context dependent.

Responses in the anterior insula were driven by proximity, direction, and speed. Importantly, in the right hemisphere, anterior insula responses also exhibited an interaction between proximity and direction, such that the approach versus retreat contrast was enhanced when the circles were “near” compared with far. The anterior insula supports subjective awareness of bodily states (Craig, 2002, 2009) and is consistently engaged during threat-related processing (Nitschke et al., 2006; Simmons et al., 2006). In particular, the anterior insula is implicated in tracking threat proximity and direction during aversive anticipation (Mobbs et al., 2010; Somerville et al., 2010). Our results replicated these findings while extending them by showing that the effects of threat proximity and direction are not independent but jointly contribute to responses in the anterior insula.

In this study, we observed a proximity by direction interaction in the right amygdala, but in the opposite direction to that seen in the anterior insula: When far, direction had a weak or no effect on responses, but when near responses were greater for retreat relative to approach. In fact, the differential response to retreat versus approach became more pronounced as the circles approached each other, with approach responses decreasing with increased proximity. In paradigms investigating the independent effects of proximity and direction on threat anticipation, Somerville et al. (2010) suggested a limited role of the amygdala in tracking threat proximity, whereas Mobbs et al. (2010) observed amygdala responses that responded to the proximity and direction of threat. In a study involving virtual predators, Mobbs et al. (2007) reported increased activation in the dorsal amygdala when threat was near, whereas responses were “stronger” in the inferior-lateral amygdala with distant threats. Thus, our results more closely resemble the latter amygdala subregion. It should be noted that in previous studies, similar to the pattern of responses observed in the current study, we and others have observed amygdala deactivations during short and long periods of sustained threat (relative to safe conditions; Grupe, Wielgosz, Davidson, & Nitschke, 2016; McMenamin et al., 2014; Choi, Padmala, & Pessoa, 2012); see also (Wager et al., 2009; Pruessner et al., 2008) in case of social stress/threat.

The role of the BST in threat processing has gained increased attention in the past two decades (Shackman & Fox, 2016; Davis & Whalen, 2001), especially during conditions involving temporally extended and less predictable threats. Given the small size of the structure and its anatomical location, studying the BST with fMRI is particularly challenging. Recently, anatomical masks for both regular and higher field scanning have been published (Theiss et al., 2017; Torrisi et al., 2015; Avery et al., 2014), which should enhance the reproducibility of published findings. We analyzed BST data using an anatomical mask and unsmoothed data, which is important because nearly all studies have used some voxelwise spatial smoothing, which blends BST signals with those of adjacent territories (beyond the inherent point spread function of imaging itself), but note that smoothing within the BST was accomplished by averaging unsmoothed time series of voxels within the anatomically defined ROI. In the right BST, we observed an effect of proximity and a proximity by speed interaction (but note that these effects were less robust as they would not survive correction for the seven tests used or 14 if one were to consider both hemispheres). The observed effect of proximity is consistent with previous findings that the BST responds to threat proximity (independent of direction; Somerville et al., 2010; Mobbs et al., 2010; although the activated region was sufficiently large as to make anatomical localization challenging in the study by Mobbs and colleagues).

The PAG of the midbrain has been implicated in aversive and defensive reactions (Bandler & Shipley, 1994; Bandler, 1988), in line with more recent studies (Tovote et al., 2016). In humans, the PAG has been suggested to be involved in negative emotional processing more generally (Satpute et al., 2013; Lindquist, Wager, Kober, Bliss-Moreau, & Barrett, 2012). The virtual tarantula manipulation by Mobbs et al. (2010), where participants were shown a prerecorded video of a spider moving toward or away from their feet was particularly effective in engaging the PAG when threat was proximal (although the activation was very extensive and thus difficult to localize). Here, in the voxelwise analysis, we only detected an effect of direction in the right midbrain/PAG where stronger responses were observed when circles were approaching compared with retreating. However, exploratory analyses revealed a robust proximity by direction interaction, as well as a proximity by direction by speed interaction. These results are potentially important because they suggest that threat-related responses in the PAG are sensitive to multiple factors that jointly determine the PAG's activity. Interestingly, unlike in the amygdala and anterior insula where we only observed an interaction between proximity and direction, speed also played a role in the PAG. However, given the exploratory nature of our analysis, future converging findings are needed to more precisely delineate the role of multiple threat-related factors on PAG activity during aversive anticipation.

A limitation of this study was that it did not include two types of control condition. First, only aversive events were encountered and not motivationally positive ones. Thus, the extent to which signals investigated here were linked to threat and not “motivational significance” more generally needs to be further investigated. Second, because a “no-shock condition” was not included, it is possible that signal fluctuations were due to processes linked to tracking circle movement, including predicting future circle positions based on current position and prior movement statistics. In this context, the anterior insula is an interesting case because it is a highly functionally diverse region and is sensitive to a very broad range of influences (Anderson, Kinnison, & Pessoa, 2013). But because anterior insula signals were sensitive to interactions between proximity and direction, it is unlikely that prediction/updating processes explained responses, as participants presumably engaged in such processing in a similar fashion when the circles were closer or father. In addition, we observed a positive correlation between proximity by direction interaction responses and responses evoked to physical shock, consistent with the fact that responses were at least in part related to anticipation of the aversive event. Finally, and more generally, the three-way interaction in the PAG (but see Results for the exploratory aspect of this result) exhibited a degree of specificity (compare left and right panels in Figure 11B) that are difficult to explain by visuocognitive processes of circle movement tracking.

Another limitation of the preset study was that participants did not have control over the threat. Unlike active avoidance paradigms where participants could perform instrumental actions to terminate or completely avoid the threat (for instance, see Mobbs et al., 2007), the passive nature of our task likely constrained the types of “defensive processing” observed. In particular, investigation of a richer set of behaviors and brain responses, such as described in the threat imminence continuum framework (Fanselow & Lester, 1988), will require novel approaches and experimental designs attuned to findings in ethology and behavioral ecology (see Mobbs, Trimmer, Blumstein, & Dayan, 2018). Finally, our choice of using nonpainful aversive stimulation was motivated by our goal to minimize potential harm to participants, and using painful stimulation likely would have generated stronger threat-related responses. Of note, a recent meta-analysis reported a large number of shared neural substrates during the processing of nonpainful and painful aversive stimuli (Hayes & Northoff, 2012).

To conclude, we investigated how multiple threat-related factors (proximity, direction, and speed) interact when varied continuously. In particular, we asked whether signal fluctuations in brain regions track threat-related factors dynamically? If so, to what factor(s) and factor combinations are they sensitive? We observed a proximity by direction interaction in the anterior insula where approach versus retreat responses were enhanced when threat was proximal. In the right amygdala, we also observed a proximity by direction interaction, but in the opposite direction as that found for the anterior insula; retreat responses were stronger than approach responses when threat was proximal. In the right BST, we observed an effect of proximity and in the right PAG/midbrain we observed an effect of direction as well as a proximity by direction by speed interaction (the latter was detected in exploratory analyses but not in a voxelwise fashion). Overall, this study refines our understanding of the mechanisms involved during aversive anticipation in the typical human brain. Importantly, it emphasizes that threat processing should be understood in a manner that is both context sensitive and dynamic. As aberrations in aversive anticipation are believed to play a major role in disorders such as anxiety and depression (Dillon et al., 2014; Grupe & Nitschke, 2013), our findings of interactions between multiple threat-related factors in regions such as the amygdala, anterior insula, and PAG may inform the understanding of brain mechanisms that are dysregulated in these disorders.

## Acknowledgments

We would like to thank Brenton McMenamin for paradigm development, Dan Levitas for data collection, Jason Smith for help with processing scripts, Mahshid Najafi for help with initial pre-processing of the fMRI data, and Nicole Friedman and Jessica Berman for help with participant recruitment. The authors also acknowledge the Behavioral and Social Sciences College, University of Maryland, high-performance computing resources (http://bsos.umd.edu/oacs/bsos-high-performance) made available for conducting the research reported in this article. The authors acknowledge funding from the National Institute of Mental Health (R01 MH071589 and R01 MH112517) and a National Science Foundation Graduate Research Fellowship to S. P.

Reprint requests should be sent to Luiz Pessoa, Department of Psychology, University of Maryland, 1147 Biology-Psychology Building, College Park, MD 20742, or via e-mail: pessoa@umd.edu.

## REFERENCES

REFERENCES
Anderson
,
M. L.
,
Kinnison
,
J.
, &
Pessoa
,
L.
(
2013
).
Describing functional diversity of brain regions and brain networks
.
Neuroimage
,
73
,
50
58
.
Avants
,
B. B.
,
Tustison
,
N. J.
, &
Song
,
G.
(
2009
).
.
Insight Journal
,
2
,
1
35
.
Avery
,
S. N.
,
Clauss
,
J. A.
,
Winder
,
D. G.
,
Woodward
,
N.
,
Heckers
,
S.
, &
Blackford
,
J. U.
(
2014
).
BNST neurocircuitry in humans
.
Neuroimage
,
91
,
311
323
.
Bach
,
D. R.
,
Flandin
,
G.
,
Friston
,
K. J.
, &
Dolan
,
R. J.
(
2009
).
Time-series analysis for rapid event-related skin conductance responses
.
Journal of Neuroscience Methods
,
184
,
224
234
.
Bandler
,
R.
(
1988
).
Brain mechanisms of aggression as revealed by electrical and chemical stimulation: Suggestion of a central role for the midbrain periaqueductal grey region
. In
A. N.
Epstein
&
A. R.
Morrison
(Eds.),
Progress in psychobiology and physiological psychology
(
Vol. 13
, pp.
67
154
).
San Diego
:
.
Bandler
,
R.
, &
Shipley
,
M. T.
(
1994
).
Columnar organization in the midbrain periaqueductal gray: Modules for emotional expression?
Trends in Neurosciences
,
17
,
379
389
.
Benjamin
,
D. J.
,
Berger
,
J. O.
,
Johannesson
,
M.
,
Nosek
,
B. A.
,
Wagenmakers
,
E.-J.
,
Berk
,
R.
, et al
(
2018
).
Redefine statistical significance
.
Nature Human Behaviour
,
2
,
6
10
.
Blanchard
,
R. J.
, &
Blanchard
,
D. C.
(
1990
).
Anti-predator defense as models of animal fear and anxiety
. In
P. F.
Brain
,
S.
Parmigiani
,
R. J.
Blanchard
, &
D.
Mainardi
(Eds.),
Ettore Majorana international life sciences series, Vol. 8. Fear and defence
(pp.
89
108
).
Amsterdam
:
.
Blanchard
,
D. C.
,
Griebel
,
G.
,
Pobbe
,
R.
, &
Blanchard
,
R. J.
(
2011
).
Risk assessment as an evolved threat detection and analysis process
.
Neuroscience & Biobehavioral Reviews
,
35
,
991
998
.
Böcker
,
K. B. E.
,
Baas
,
J. M. P.
,
Kenemans
,
J. L.
, &
Verbaten
,
M. N.
(
2001
).
Stimulus-preceding negativity induced by fear: A manifestation of affective anticipation
.
International Journal of Psychophysiology
,
43
,
77
90
.
Brown
,
C. A.
,
Seymour
,
B.
,
Boyle
,
Y.
,
El-Deredy
,
W.
, &
Jones
,
A. K. P.
(
2008
).
Modulation of pain ratings by expectation and uncertainty: Behavioral characteristics and anticipatory neural correlates
.
Pain
,
135
,
240
250
.
Choi
,
J. M.
,
,
S.
, &
Pessoa
,
L.
(
2012
).
Impact of state anxiety on the interaction between threat monitoring and cognition
.
Neuroimage
,
59
,
1912
1923
.
Cohen
,
M. S.
(
1997
).
Parametric analysis of fMRI data using linear systems methods
.
Neuroimage
,
6
,
93
103
.
Cox
,
R. W.
(
1996
).
AFNI: Software for analysis and visualization of functional magnetic resonance neuroimages
.
Computers and Biomedical Research
,
29
,
162
173
.
Cox
,
R. W.
,
Chen
,
G.
,
Glen
,
D. R.
,
Reynolds
,
R. C.
, &
Taylor
,
P. A.
(
2017
).
fMRI clustering in AFNI: False-positive rates redux
.
Brain Connectivity
,
7
,
152
171
.
Craig
,
A. D.
(
2002
).
How do you feel? Interoception: The sense of the physiological condition of the body
.
Nature Reviews Neuroscience
,
3
,
655
666
.
Craig
,
A. D.
(
2009
).
How do you feel—Now? The anterior insula and human awareness
.
Nature Reviews Neuroscience
,
10
,
59
70
.
Davis
,
M.
,
Walker
,
D. L.
,
Miles
,
L.
, &
Grillon
,
C.
(
2010
).
Phasic vs. sustained fear in rats and humans: Role of the extended amygdala in fear vs. anxiety
.
Neuropsychopharmacology
,
35
,
105
135
.
Davis
,
M.
, &
Whalen
,
P. J.
(
2001
).
The amygdala: Vigilance and emotion
.
Molecular Psychiatry
,
6
,
13
34
.
Dillon
,
D. G.
,
Rosso
,
I. M.
,
Pechtel
,
P.
,
Killgore
,
W. D. S.
,
Rauch
,
S. L.
, &
Pizzagalli
,
D. A.
(
2014
).
Peril and pleasure: An RDoC-inspired examination of threat responses and reward processing in anxiety and depression
.
Depression & Anxiety
,
31
,
233
249
.
Eklund
,
A.
,
Nichols
,
T. E.
, &
Knutsson
,
H.
(
2016
).
Cluster failure: Why fMRI inferences for spatial extent have inflated false-positive rates
.
Proceedings of the National Academy of Sciences, U.S.A.
,
113
,
7900
7905
.
Engelmann
,
J. B.
,
Meyer
,
F.
,
Fehr
,
E.
, &
Ruff
,
C. C.
(
2015
).
Anticipatory anxiety disrupts neural valuation during risky choice
.
Journal of Neuroscience
,
35
,
3085
3099
.
Fanselow
,
M. S.
, &
Lester
,
L. S.
(
1988
).
A functional behavioristic approach to aversively motivated behavior: Predatory imminence as a determinant of the topography of defensive behavior
. In
R. C.
Bolles
&
M. D.
Beecher
(Eds.),
Evolution and learning
(pp.
185
212
).
Hillsdale, NJ
:
Erlbaum
.
Feinberg
,
D. A.
,
Moeller
,
S.
,
Smith
,
S. M.
,
Auerbach
,
E.
,
Ramanna
,
S.
,
Glasser
,
M. F.
, et al
(
2010
).
Multiplexed echo planar imaging for sub-second whole brain fMRI and fast diffusion imaging
.
PLoS One
,
5
,
e15710
.
Fox
,
A. S.
,
Oler
,
J. A.
,
Tromp
,
D. P. M.
,
Fudge
,
J. L.
, &
Kalin
,
N. H.
(
2015
).
Extending the amygdala in theories of threat processing
.
Trends in Neurosciences
,
38
,
319
329
.
Greve
,
D. N.
, &
Fischl
,
B.
(
2009
).
Accurate and robust brain image alignment using boundary-based registration
.
Neuroimage
,
48
,
63
72
.
Grupe
,
D. W.
, &
Nitschke
,
J. B.
(
2013
).
Uncertainty and anticipation in anxiety: An integrated neurobiological and psychological perspective
.
Nature Reviews Neuroscience
,
14
,
488
501
.
Grupe
,
D. W.
,
Wielgosz
,
J.
,
Davidson
,
R. J.
, &
Nitschke
,
J. B.
(
2016
).
Neurobiological correlates of distinct post-traumatic stress disorder symptom profiles during threat anticipation in combat veterans
.
Psychological Medicine
,
46
,
1885
1895
.
Hayes
,
D. J.
, &
Northoff
,
G.
(
2012
).
Common brain activations for painful and non-painful aversive stimuli
.
BMC Neuroscience
,
13
,
60
.
Iglesias
,
J. E.
,
Liu
,
C.-Y.
,
Thompson
,
P. M.
, &
Tu
,
Z.
(
2011
).
Robust brain extraction across datasets and comparison with publicly available methods
.
IEEE Transactions on Medical Imaging
,
30
,
1617
1634
.
Jensen
,
J.
,
McIntosh
,
A. R.
,
Crawley
,
A. P.
,
Mikulis
,
D. J.
,
Remington
,
G.
, &
Kapur
,
S.
(
2003
).
Direct activation of the ventral striatum in anticipation of aversive stimuli
.
Neuron
,
40
,
1251
1257
.
Kavaliers
,
M.
, &
Choleris
,
E.
(
2001
).
Antipredator responses and defensive behavior: Ecological and ethological approaches for the neurosciences
.
Neuroscience & Biobehavioral Reviews
,
25
,
577
586
.
Lim
,
C. L.
,
Rennie
,
C.
,
Barry
,
R. J.
,
Bahramali
,
H.
,
Lazzaro
,
I.
,
Manor
,
B.
, et al
(
1997
).
Decomposing skin conductance into tonic and phasic components
.
International Journal of Psychophysiology
,
25
,
97
109
.
Lindquist
,
K. A.
,
Wager
,
T. D.
,
Kober
,
H.
,
Bliss-Moreau
,
E.
, &
Barrett
,
L. F.
(
2012
).
The brain basis of emotion: A meta-analytic review
.
Behavioral and Brain Sciences
,
35
,
121
143
.
McMenamin
,
B. W.
,
Langeslag
,
S. J. E.
,
Sirbu
,
M.
,
,
S.
, &
Pessoa
,
L.
(
2014
).
Network organization unfolds over time during periods of anxious anticipation
.
Journal of Neuroscience
,
34
,
11261
11273
.
McShane
,
B. B.
,
Gal
,
D.
,
Gelman
,
A.
,
Robert
,
C.
, &
Tackett
,
J. L.
(
2017
).
Abandon statistical significance
.
arXiv preprint arXiv:170907588
.
Mobbs
,
D.
,
Hagan
,
C. C.
,
Dalgleish
,
T.
,
Silston
,
B.
, &
Prévost
,
C.
(
2015
).
The ecology of human fear: Survival optimization and the nervous system
.
Frontiers in Neuroscience
,
9
,
55
.
Mobbs
,
D.
,
Petrovic
,
P.
,
Marchant
,
J. L.
,
Hassabis
,
D.
,
Weiskopf
,
N.
,
Seymour
,
B.
, et al
(
2007
).
When fear is near: Threat imminence elicits prefrontal-periaqueductal gray shifts in humans
.
Science
,
317
,
1079
1083
.
Mobbs
,
D.
,
Trimmer
,
P. C.
,
Blumstein
,
D. T.
, &
Dayan
,
P.
(
2018
).
Foraging for foundations in decision neuroscience: Insights from ethology
.
Nature Reviews Neuroscience
,
19
,
419
427
.
Mobbs
,
D.
,
Yu
,
R.
,
Rowe
,
J. B.
,
Eich
,
H.
,
FeldmanHall
,
O.
, &
Dalgleish
,
T.
(
2010
).
Neural activity associated with monitoring the oscillating threat value of a tarantula
.
Proceedings of the National Academy of Sciences, U.S.A.
,
107
,
20582
20586
.
Mumford
,
J. A.
,
Poline
,
J.-B.
, &
Poldrack
,
R. A.
(
2015
).
Orthogonalization of regressors in fMRI models
.
PLoS One
,
10
,
e0126255
.
Najafi
,
M.
,
Kinnison
,
J.
, &
Pessoa
,
L.
(
2017
).
Dynamics of intersubject brain networks during anxious anticipation
.
Frontiers in Human Neuroscience
,
11
,
552
.
Neter
,
J.
,
Kutner
,
M. H.
,
Nachtsheim
,
C. J.
, &
Wasserman
,
W.
(
1996
).
Applied linear statistical models
(4th ed.).
Chicago
:
Irwin
.
Nitschke
,
J. B.
,
Sarinopoulos
,
I.
,
Mackiewicz
,
K. L.
,
Schaefer
,
H. S.
, &
Davidson
,
R. J.
(
2006
).
Functional neuroanatomy of aversion and its anticipation
.
Neuroimage
,
29
,
106
116
.
Pessoa
,
L.
(
2016
).
The emotional brain
. In
P. M.
Conn
(Ed.),
Conn's translational neuroscience
(pp.
635
656
).
Amsterdam
:
Elsevier
.
Pruessner
,
J. C.
,
Dedovic
,
K.
,
Khalili-Mahani
,
N.
,
Engert
,
V.
,
Pruessner
,
M.
,
Buss
,
C.
, et al
(
2008
).
Deactivation of the limbic system during acute psychosocial stress: Evidence from positron emission tomography and functional magnetic resonance imaging studies
.
Biological Psychiatry
,
63
,
234
240
.
Rousselet
,
G. A.
, &
Pernet
,
C. R.
(
2012
).
Improving standards in brain-behavior correlation analyses
.
Frontiers in Human Neuroscience
,
6
,
119
.
Satpute
,
A. B.
,
Wager
,
T. D.
,
,
J.
,
Bianciardi
,
M.
,
Choi
,
J.-K.
,
Buhle
,
J. T.
, et al
(
2013
).
Identification of discrete functional subregions of the human periaqueductal gray
.
Proceedings of the National Academy of Sciences, U.S.A.
,
110
,
17101
17106
.
Shackman
,
A. J.
, &
Fox
,
A. S.
(
2016
).
Contributions of the central extended amygdala to fear and anxiety
.
Journal of Neuroscience
,
36
,
8050
8063
.
Shattuck
,
D. W.
, &
Leahy
,
R. M.
(
2002
).
BrainSuite: An automated cortical surface identification tool
.
Medical Image Analysis
,
6
,
129
142
.
Siegel
,
J. S.
,
Power
,
J. D.
,
Dubis
,
J. W.
,
Vogel
,
A. C.
,
Church
,
J. A.
,
Schlaggar
,
B. L.
, et al
(
2014
).
Statistical improvements in functional magnetic resonance imaging analyses produced by censoring high-motion data points
.
Human Brain Mapping
,
35
,
1981
1996
.
Simmons
,
A.
,
Strigo
,
I.
,
Matthews
,
S. C.
,
Paulus
,
M. P.
, &
Stein
,
M. B.
(
2006
).
Anticipation of aversive visual stimuli is associated with increased insula activation in anxiety-prone subjects
.
Biological Psychiatry
,
60
,
402
409
.
Smith
,
S. M.
,
Jenkinson
,
M.
,
Woolrich
,
M. W.
,
Beckmann
,
C. F.
,
Behrens
,
T. E. J.
,
Johansen-Berg
,
H.
, et al
(
2004
).
Advances in functional and structural MR image analysis and implementation as FSL
.
Neuroimage
,
23(Suppl. 1)
,
S208
S219
.
Somerville
,
L. H.
,
Whalen
,
P. J.
, &
Kelley
,
W. M.
(
2010
).
Human bed nucleus of the stria terminalis indexes hypervigilant threat monitoring
.
Biological Psychiatry
,
68
,
416
424
.
Spielberger
,
C. D.
,
Gorsuch
,
R. L.
, &
Lushene
,
R. E.
(
1970
).
Manual for the State–Trait Anxiety Inventory
.
Palo Alto, CA
:
Consulting Psychologists Press
.
Theiss
,
J. D.
,
Ridgewell
,
C.
,
McHugo
,
M.
,
Heckers
,
S.
, &
Blackford
,
J. U.
(
2017
).
Manual segmentation of the human bed nucleus of the stria terminalis using 3 T MRI
.
Neuroimage
,
146
,
288
292
.
Torrisi
,
S.
,
O'Connell
,
K.
,
Davis
,
A.
,
Reynolds
,
R.
,
Balderston
,
N.
,
Fudge
,
J. L.
, et al
(
2015
).
Resting state connectivity of the bed nucleus of the stria terminalis at ultra-high field
.
Human Brain Mapping
,
36
,
4076
4088
.
Tovote
,
P.
,
Esposito
,
M. S.
,
Botta
,
P.
,
Chaudun
,
F.
,
,
J. P.
,
Markovic
,
M.
, et al
(
2016
).
Midbrain circuits for defensive behaviour
.
Nature
,
534
,
206
212
.
Vytal
,
K. E.
,
Overstreet
,
C.
,
Charney
,
D. R.
,
Robinson
,
O. J.
, &
Grillon
,
C.
(
2014
).
Sustained anxiety increases amygdala–dorsomedial prefrontal coupling: A mechanism for maintaining an anxious state in healthy adults
.
Journal of Psychiatry & Neuroscience
,
39
,
321
329
.
Wager
,
T. D.
,
Waugh
,
C. E.
,
Lindquist
,
M.
,
Noll
,
D. C.
,
Fredrickson
,
B. L.
, &
Taylor
,
S. F.
(
2009
).
Brain mediators of cardiovascular responses to social threat: Part I. Reciprocal dorsal and ventral sub-regions of the medial prefrontal cortex and heart-rate reactivity
.
Neuroimage
,
47
,
821
835
.
Wilcox
,
R. R.
(
2012
).
Introduction to robust estimation and hypothesis testing
(3rd ed.).
Cambridge, MA
:
.
Woo
,
C.-W.
,
Krishnan
,
A.
, &
Wager
,
T. D.
(
2014
).
Cluster-extent based thresholding in fMRI analyses: Pitfalls and recommendations
.
Neuroimage
,
91
,
412
419
.

## Author notes

*

These authors contributed equally to this work.