A Shared Threat-Anticipation Circuit Is Dynamically Engaged at Different Moments by Certain and Uncertain Threat

Paradigm structure and design considerations

The Maryland Threat Countdown paradigm is a well-established, fMRI-optimized variant of temporally uncertain-threat assays that have been validated using fear-potentiated startle and acute anxiolytic administration (e.g., benzodiazepine) in mice, rats, and humans (Miles et al., 2011; Hefner et al., 2013; Daldrup et al., 2015; Lange et al., 2017; Moberg et al., 2017). The paradigm has been successfully deployed in independent samples of samples of university students and community volunteers (Kim et al., 2023; Grogans et al., 2024).

As shown schematically in Figure 1, the paradigm takes the form of a 2 (Valence: Threat, Safety) × 2 (Temporal Certainty: Uncertain, Certain) randomized, event-related, repeated-measures design (3 scans; 6 trials/condition/scan). Participants were completely informed about the task design and contingencies prior to scanning. Simulations were used to optimize the detection and deconvolution of task-related hemodynamic signals. Stimulus presentation was controlled using Presentation software (version 19.0, Neurobehavioral Systems).

On certain-threat trials, participants saw a descending stream of integers (“countdown;” e.g., 30, 29, 28...3, 2, 1) for 18.75 s. To ensure robust fear and anxiety, this anticipation epoch culminated with the presentation of a noxious electric shock, unpleasant photograph (e.g., mutilated body), and thematically related audio clip (e.g., scream). Uncertain-threat trials were similar, but the integer stream was randomized and presented for an uncertain and variable duration (8.75–30.00 s; M = 18.75 s). Participants knew that something aversive was going to occur, but they had no way of knowing precisely when. Consistent with methodological recommendations (Shackman and Fox, 2016), the average duration of the anticipation epoch was identical across conditions, ensuring an equal number of measurements (TRs/condition). The specific mean duration was chosen to enhance detection of task-related differences in the blood oxygen level-dependent (BOLD) signal (“activation”; Henson, 2007) and to allow sufficient time for sustained responses to become evident. Certain- and uncertain-safety trials were similar but terminated with the presentation of benign reinforcers (see below). Valence was continuously signaled by the background color of the display. Temporal certainty was signaled by the nature of the integer stream. Certain trials always began with the presentation of the number 30. On uncertain trials, integers were randomly drawn from a near-uniform distribution ranging from 1 to 45 to reinforce the impression that they could be much shorter or longer than certain trials and to minimize incidental temporal learning (“time-keeping”). To demonstrate the variable duration of Uncertain trials, during scanning, the first three uncertain trials featured short (8.75 s), medium (15.00 s), and long (28.75 s) anticipation epochs. To mitigate potential confusion and eliminate mnemonic demands, a lower-case “c” or “u” was presented at the lower edge of the display throughout the anticipatory epoch. White-noise visual masks (3.2 s) were presented between trials to minimize the persistence of visual reinforcers in iconic memory. Anticipatory distress ratings and skin conductance were also acquired, as detailed previously (Grogans et al., 2024).

Procedures

Prior to scanning, participants practiced an abbreviated version of the paradigm (without electrical stimulation) until they indicated and staff confirmed understanding. Benign and aversive electrical stimulation levels were individually titrated.

Electrical stimuli

Electrical stimuli (100 ms; 2 ms pulses every 10 ms) were generated using an MRI-compatible constant-voltage stimulator system (STMEPM-MRI; Biopac Systems) and delivered using MRI-compatible, disposable carbon electrodes (Biopac) attached to the fourth and fifth digits of the nondominant hand.

Visual stimuli

A total of 72 aversive and benign photographs (1.8 s) were selected from the International Affective Picture System (Hur et al., 2020b). Visual stimuli were digitally back-projected (Powerlite Pro G5550, Epson America) onto a semi-opaque screen mounted at the head-end of the scanner bore and viewed using a mirror mounted on the head-coil.

Auditory stimuli

Seventy-two aversive and benign auditory stimuli (0.8 s) were adapted from open-access online sources and delivered using an amplifier (PA-1 Whirlwind) with in-line noise-reducing filters and ear buds (S14; Sensimetrics) fitted with noise-reducing ear plugs (Hearing Components).

Anatomical data processing

T1- and T2-weighted images were inhomogeneity corrected using N4 (Tustison et al., 2010) and denoised using ANTS (Avants et al., 2011). The brain was then extracted using BEaST (Eskildsen et al., 2012) and brain-extracted and normalized reference brains from IXI (BIAC, 2022). Brain-extracted T1 images were normalized to a version of the brain-extracted 1 mm T1-weighted MNI152 template (nonlinear 6th-generation symmetric average; Grabner et al., 2006) modified to remove extracerebral tissue. Normalization was performed using the diffeomorphic approach implemented in SyN (version 2.3.4; Avants et al., 2011). T2-weighted images were rigidly coregistered with the corresponding T1 prior to normalization. The brain extraction mask from the T1 was then applied. Tissue priors were unwarped to native space using the inverse of the diffeomorphic transformation (Lorio et al., 2016). Brain-extracted T1 and T2 images were segmented using native-space priors generated in FAST (version 6.0.4; Jenkinson et al., 2012) for subsequent use in T1-EPI coregistration (see below).

Fieldmap data processing

SE images and topup were used to create fieldmaps. Fieldmaps were converted to radians, median-filtered, and smoothed (2 mm). The average of the distortion-corrected SE images was inhomogeneity corrected using N4 and masked to remove extracerebral voxels using 3dSkullStrip (AFNI version 23.1.10). The resulting mask was minimally eroded to further exclude extracerebral voxels.

Functional data processing

EPI files were despiked (3dDespike), slice time corrected to the TR-center using 3dTshift, and motion-corrected to the first volume and inhomogeneity corrected using ANTS (12-parameter affine). Transformations were saved in ITK-compatible format for subsequent processing (McCormick et al., 2014). The first volume was extracted for EPI-T1 coregistration. The reference EPI volume was simultaneously coregistered with the corresponding T1-weighted image in native space and corrected for geometric distortions using boundary-based registration (Jenkinson et al., 2012). This step incorporated the previously created fieldmap, undistorted SE, T1, white matter (WM) image, and masks. The spatial transformations necessary to transform each EPI volume from native space to the reference EPI, from the reference EPI to the T1, and from the T1 to the template were concatenated and applied to the processed EPI data in a single step to minimize incidental spatial blurring. Normalized EPI data were resampled (2 mm³) using fifth-order b-splines. Voxelwise analyses employed data that were spatially smoothed (4 mm) using 3DblurInMask. To minimize signal mixing, smoothing was confined to the gray-matter compartment, defined using a variant of the Harvard–Oxford cortical and subcortical atlases that was expanded to include the bed nucleus of the stria terminalis (BST) and periaqueductal gray (PAG; Frazier et al., 2005; Desikan et al., 2006; Makris et al., 2006; Edlow et al., 2012; Theiss et al., 2017). Focal analyses of the EAc leveraged spatially unsmoothed data and anatomically defined regions of interest (see below), as in prior work (Tillman et al., 2018; Kim et al., 2023; Grogans et al., 2024; Hur et al., 2025).

Data exclusions

Volume-to-volume displacement (>0.5 mm) was used to assess residual motion artifact. Scans with excessively frequent residual artifacts (>2 SD) were discarded. Participants with insufficient usable fMRI data (<2 scans) were excluded from analyses (see above).

Overview of first-level (single-subject) fMRI modeling

For each participant, first-level modeling was performed using general linear models (GLMs) implemented in 3dREMLfit (ARMA_1,1; fourth-order Legendre high-pass filter). Regressors were convolved with the SPM12 canonical hemodynamic-response function (HRF). Epochs corresponding to the presentation of the four types of reinforcers, white-noise visual masks, and rating prompts were simultaneously modeled using the same approach. As in our prior work, nuisance variates included volume-to-volume displacement and its first derivative, six motion parameters and their first derivatives, cerebrospinal fluid (CSF) signal, instantaneous pulse and respiration rates, and nuisance signals (e.g., brain edge, CSF edge, global motion, white matter, extracerebral soft tissue; Anderson et al., 2011; Pruim et al., 2015). Volumes with excessive volume-to-volume displacement (>0.75 mm) and those during and immediately following reinforcer delivery were censored. EPI volumes acquired before the first trial and following the final trial were unmodeled and contributed to the baseline estimate.

Conventional “boxcar” model

The present sample of 220 datasets represents a superset of the 99 featured in an earlier report from our group that employed a conventional “boxcar” fMRI modeling approach and an older data-processing pipeline (Hur et al., 2020b). As a precursor to hypothesis testing, we used a conventional first-level model to confirm that the larger, reprocessed dataset broadly reproduced our published observations. Hemodynamic reactivity to the threat-anticipation paradigm was modeled using variable-duration rectangular (“boxcar”) regressors that spanned the entirety of the anticipation (“countdown”) epoch for uncertain-threat, certain-threat, and uncertain-safety trials (8.75–30.00 s; Fig. 1). To maximize design efficiency, certain-safety anticipation served as the first-level reference condition and contributed to the baseline estimate (Poline et al., 2007).

Onset–Sustained–Phasic model

Neuroimaging research by our group and others has relied on simplified “boxcar” modeling approaches that reduce the neural dynamics anticipated by theory and psychophysiological research to a single average response (Wang et al., 2024). Here we used two complementary hemodynamic models to characterize the time-varying signals elicited by certain- and uncertain-threat anticipation. The Onset–Sustained–Phasic (OSP) model used a multiple-regression framework to identify the variance in threat-anticipation signals that was uniquely associated, in the partial-correlation sense, with temporally overlapping Onset, Sustained, and Phasic regressors (Fig. 2a). The first-level design matrix incorporated a punctate event or “impulse” time locked to the onset of the anticipation epoch, a variable-duration rectangular function that spanned the entirety of the anticipation epoch (to capture sustained increases in activation), and a rectangular-function time locked to the offset of the anticipation epoch (to capture phasic surges in activation just prior to threat encounters). Paralleling the “boxcar” model, the Sustained regressor for certain-safety served as the first-level reference condition and contributed to the baseline estimate. The duration of the Phasic regressor (6.15 s) was chosen based on a combination of theory and simulations aimed at minimizing regressor colinearity. The mean condition-wise variance inflation factor was <1.93 for the full task-related model (excluding nuisance regressors). As detailed in Figure 2, the Sustained regressor captures variance in the hemodynamic signal associated with a particular trial type (e.g., uncertain-threat anticipation) above-and-beyond that captured by the Onset and Phasic regressors. Unlike conventional boxcar models, this provides an estimate of sustained activation that is unconfounded by nonspecific orienting or salience responses reflexively triggered by the dramatic change in sensory stimulation associated with trial onset (Sokolov et al., 2002; Menon, 2015). Likewise, the Phasic regressor captures variance in the hemodynamic signal above-and-beyond that captured by the Onset and Sustained regressors, with positive coefficients indicating an increase in activation in the final moments of the anticipation epoch relative to that associated with the Sustained and Onset regressors. In sum, the OSP model casts the overall magnitude of the hemodynamic signal as a linear combination of the Onset, Sustained, and Phasic regressors; nuisance regressors; and error (Fig. 2).

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

Theory-driven hemodynamic modeling. We used two modeling approaches to investigate time-varying responses to certain- and uncertain-threat anticipation. The relative merits of the two models are described in the main text (Materials and Methods, Complementary strengths and limitations of the two models). Because these models are uncommon, we provide a general overview in panels a and b and, for readers interested in greater detail, a description of their nuances in panels c–f. a, Overview of the OSP model. The OSP model used multiple regression to partition the variance in threat-anticipation signals uniquely associated with temporally overlapping Onset, Sustained, and Phasic regressors (see main text for details). The OSP is ideal for examining between within-moment contrasts (e.g., Phasic regressor: certain vs uncertain threat). b, Overview of the Convolved Blocks (CB) model. The CB model splits the anticipation epoch into a sequence of short (6.25 s), nonoverlapping rectangular functions or “blocks,” each convolved with a canonical HRF. The CB model is ideal for examining between-moment contrasts and temporal trends (e.g., quadratic effects), and it enabled a rigorous test of surges in activation in the moments just prior to certain-threat encounters (Certain-Threat: late vs middle block)—something not permitted by the OSP model. c, Consider a voxel that shows a sustained level of heightened activation during certain-threat anticipation (red line). In the OSP model, this is captured by a strong loading or “weight” on the Sustained regressor and nil loadings on the Onset and Phasic regressors. d, In the CB model, this pattern of weights is instead associated with a transient increase in the hemodynamic signal in the middle of the anticipation epoch. f, In the CB model, sustained activation is captured by uniformly strong loadings across the early, middle, and late regressors. e, In the OSP model, this pattern of weights is instead associated with transient onset and phasic responses, superimposed on a strong sustained response. In effect, the Onset and Phasic regressors serve to modulate the leading and trailing edges of a sustained wave of activation. It merits comment that, although the three weights are equally strong (barplot: O ≈ S ≈ P), the moment-by-moment height of the hemodynamic signal is not (red line: O > S < P). This reflects the fact that the OSP model casts the threat-anticipation signal as the linear combination of 3 temporally overlapping regressors. g, In the OSP model, the relative height of activation (red line: S < P) can be reversed from the rank order of the weights (barplot: S > P). h, In contrast, the CB model provides a one-to-one mapping between the moment-by-moment height of the hemodynamic signal and the early, middle, and late regression weights. BOLD, blood oxygenation level-dependent; CB, Convolved Blocks model; HRF, hemodynamic response function; fMRI, functional magnetic resonance imaging; OSP, onset-sustained-phasic; s, seconds.

Convolved Blocks model

To clarify interpretation (see next section), we employed a piecewise approach that arbitrarily splits the anticipation epoch into a sequence of 2–5 short (6.25 s), nonoverlapping rectangular functions or “blocks”, each convolved with a canonical HRF (Fig. 2b). Here the second certain-safety block (6.25–12.5 s) served as the first-level reference condition and contributed to the baseline estimate. Regressor colinearity was acceptable (Mumford et al., 2015). The mean condition-wise variance inflation factor was <1.91 for the full task-related model.

Complementary strengths and limitations of the two models

The OSP and Convolved Blocks models have complementary strengths and limitations (Fig. 2). The OSP Phasic regressor is time locked to the offset of the anticipation epoch, ensuring that it always indexes neural activation in the seconds just before threat is encountered—regardless of its temporal certainty or the overall duration of the anticipation (“countdown”) epoch. Because it is effectively a partial correlation, the OSP Phasic regressor captures variance in activation above-and-beyond that captured by the temporally overlapping Sustained regressor, providing a more “pure” or conservative estimate of phasic surges in activation (Fig. 2c,e,g). In short, the OSP model provides a unified, nonarbitrary way to address the variable duration of different trial types, making it ideal for examining within-moment contrasts (e.g., Phasic regressor: certain vs uncertain threat). Despite these strengths, the partial correlations yielded by OSP model do not permit straightforward assessments of dynamic changes in activation. In contrast, the Convolved Blocks model, while arbitrary in timing, yields activation estimates that are statistically independent and directly comparable across moments of time (Fig. 2d,f,h). As such, the Convolved Blocks model is ideal for examining between-moment contrasts and overall temporal trends in anticipatory activation (e.g., quadratic effects). In particular, the Convolved Blocks model provides a natural way to rigorously test hypothesized surges in activation in the seconds just prior to temporally certain encounters with threat (Certain-Threat: late vs middle block).

EAc ROIs

The central EAc occupies center stage in most neurobiological models of fear and anxiety, including RDoC (NIMH, 2011, 2020a,b; Tovote et al., 2015; LeDoux and Pine, 2016; Fox and Shackman, 2019; Mobbs et al., 2019; Tseng et al., 2023; Akiki et al., 2025). The EAc is a functional macrocircuit encompassing the central nucleus of the amygdala (Ce) and the neighboring BST (Fox et al., 2015b; Shackman and Fox, 2016; Fox and Shackman, 2019). As in prior work by our group (Grogans et al., 2024), EAc activation was quantified using anatomically defined probabilistic ROIs (Theiss et al., 2017; Tillman et al., 2018). The BST ROI was modified to remove voxels encroaching upon neighboring regions of the striatum, thalamus, and ventricles. It mostly encompasses the supracommissural BST, given the difficulty of reliably discriminating the borders of regions below the anterior commissure in T1-weighted images (Walter et al., 1991; Kruger et al., 2015). Bilateral ROIs were decimated to the 2 mm resolution of the fMRI data.

Participants

A total of 241 participants were recruited and scanned. Of these, 6 withdrew from the study due to excess distress during the imaging session and 1 withdrew for undisclosed reasons following the imaging session. Another 14 participants were excluded from fMRI analyses due to incidental neurological findings (n = 4), technical problems (n = 2), or insufficient usable data (n = 8; see below), yielding a racially diverse sample of 220 participants (49.5% female; 61.4% White, 18.2% Asian, 8.6% African American, 4.1% Hispanic, 7.3% Multiracial/Other; M = 18.8 years, SD = 0.4 years).

Power analysis

To enable readers to better interpret in our results, we performed a post hoc power analysis. G-Power (version 3.1.9.2) indicated that the final sample of 220 usable fMRI datasets provides 80% power to detect a benchmark (“generic”) mean difference as small as Cohen's d = 0.19 (α = 0.05, two-tailed; Cohen, 1988; Faul et al., 2007). The study was not preregistered.

Analytics overview

Analyses were performed using a combination of SPM12 (Wellcome Centre for Human Neuroimaging, 2022), SPSS (version 27.0.1), and JASP (version 0.16.4.0; Love et al., 2019). Diagnostic procedures and data visualizations were used to confirm that test assumptions were satisfied (Tukey, 1977). Some figures were created using R (version 4.0.2), RStudio (version 1.2.1335), tidyverse (version 2.0), ggplot2 (version 3.3.6), ggpubr (version 0.6.0), plotrix (version 3.8-4), and MRIcron (version 1.0.20190902; Lemon, 2006; Wickham, 2016; Rorden, 2019; Wickham et al., 2019; R Core Team, 2022; RStudio Team, 2022; Kassambara, 2023). Clusters and peaks were labeled using the Harvard–Oxford atlas, supplemented by the Mai, Allen, Jülich (version 3.1), and other atlases and resources (Frazier et al., 2005; Desikan et al., 2006; Makris et al., 2006; Mai et al., 2015; Ding et al., 2016; Theiss et al., 2017; ten Donkelaar et al., 2018; Tillman et al., 2018; Amunts et al., 2021).

Resource sharing

Raw data are available at the National Institute of Mental Health Data Archive (https://nda.nih.gov/edit_collection.html?id=2447). Neuroimaging maps are available at NeuroVault (https://neurovault.org/collections/15274). Task materials, statistical code, de-identified processed data, and neuroimaging cluster tables are available at OSF (https://osf.io/e2ngf). The negative affect brain signature is available at GitHub (https://github.com/canlab/Neuroimaging_Pattern_Masks/tree/master/Multivariate_signature_patterns/2021_Ceko_MPA2_multiaversive).

Conventional “boxcar” model

The present sample of 220 datasets represents a superset of the 99 featured in an earlier report that employed a conventional “boxcar” hemodynamic-modeling approach, older data-processing pipeline, and larger spatial-smoothing kernel (6 mm; Hur et al., 2020b). Here we used standard voxelwise GLMs to confirm that conventional modeling of the larger, reprocessed dataset broadly reproduced our published results (FDR q < 0.05, whole-brain corrected). As in prior work by our group (Hur et al., 2020b; Kim et al., 2023), a minimum-conjunction analysis (Logical “AND”) was used to identify regions showing significant activation during the anticipation of certain and uncertain threat relative to their respective second-level reference conditions (Nichols et al., 2005).

Voxelwise analyses of sustained and phasic activation dynamics

A major goal of the present study was to identify regions showing sustained levels of heighted activation during uncertain-threat anticipation and phasic surges in activation during the final moments of certain-threat anticipation. Standard voxelwise GLMs were used to compare each kind of anticipated threat (e.g., uncertain threat) with its second-level reference condition (e.g., uncertain safety) and to one another (e.g., certain threat; FDR q < 0.05, whole-brain corrected). Hypothesis testing focused on the OSP Sustained and Phasic regressors. A minimum-conjunction analysis was used to identify regions showing significant sustained activation during uncertain-threat anticipation and significant phasic activation during the terminal portion of certain-threat anticipation, that is, regions showing evidence of neuroanatomical colocalization (Fox and Shackman, 2019; Hur et al., 2020b; Shackman and Fox, 2021; Shackman et al., 2024). While not a focus of the present study, exploratory analyses of the OSP Onset regressor, which mainly captures reflexive orienting responses, are also briefly summarized.

Focused tests of phasic activation

While the OSP model is well suited for within-moment contrasts (e.g., Phasic regressor: certain vs uncertain threat), the resulting partial-regression coefficients do not permit straightforward interpretation of between-moment contrasts (for details, see above and Fig. 2). To more fully test phasic effects, we used activation estimates from the Convolved Blocks model and a standard voxelwise GLM to identify regions showing significant increases in activation during the final (12.5–18.75 s) relative to the middle (6.25–12.5 s) third of the certain-threat anticipation epoch (FDR q < 0.05, whole-brain corrected).

We used activation estimates derived from the Convolved Blocks model and a standard voxelwise GLM to identify regions that differed in their responses to the final block of certain threat (an indicator of phasic surges in activation) compared with the second block of uncertain threat (an indicator of sustained activation), each relative to their respective reference condition (FDR q < 0.05, whole-brain corrected).

EAc ROI analyses

ROI analyses used activation estimates (i.e., standardized regression coefficients) for convolved blocks 1–3, extracted and averaged for each contrast (e.g., uncertain-threat vs uncertain-safety anticipation), region, and participant. This approach enabled us to span the mean duration of the anticipation epoch (0–18.75 s) and examine overall temporal trends. Unlike conventional whole-brain voxelwise analyses—which screen thousands of voxels for statistical significance and yield optimistically biased associations—anatomically defined ROIs “fix” the measurements-of-interest a priori, providing statistically unbiased effect-size estimates (Poldrack et al., 2017). As a precursor to hypothesis testing, we used one-sample Student's t tests to confirm that the Ce and BST ROIs are nominally engaged by anticipated threat (p < 0.05, uncorrected). For hypothesis testing, we used a standard 2 (Region: BST, Ce) × 2 (Threat-Certainty: Certain, Uncertain) × 3 (Convolved Block: 1, 2, 3) repeated-measures GLM (Huynh–Feldt correction) to interrogate potential regional differences in activation dynamics. Paralleling the voxelwise analyses, hypothesis testing focused on the second block (6.25–12.5 s) of uncertain-threat anticipation (a proxy for sustained activation) and the final block (12.5–18.75 s) of certain-threat anticipation (a proxy for phasic surges in activation). Significant GLM effects were further interrogated using planned linear and quadratic polynomial contrasts, as in prior psychophysiological research (Grillon et al., 1993a,b; Löw et al., 2015). Polynomial analyses enabled us to evaluate whether the time course of activation differed as a function of Region, Threat-Certainty, or their interaction. In particular, we sought to test whether certain-and-imminent threat is associated with phasic surges in activation in the final third of the anticipation epoch (relative to the middle third), whether this hypothesized quadratic trend is stronger than that evinced during the identical moments of uncertain-threat anticipation, and whether these temporal dynamics differ between the BST and Ce.

Frequentist (Cohen's d) and Bayesian (BF₁₀) effect sizes were used to clarify nonsignificant regional differences. Cohen's d was interpreted using established benchmarks (Cohen, 1988, 1994; Schimmack, 2019), ranging from small (d = 0.20) to nil (d ≤ 0.10). BF₁₀ quantifies the relative performance of the null hypothesis (H₀; e.g., the absence of a credible mean difference) and the alternative hypothesis (H₁; e.g., the presence of a credible mean difference) on a 0 to ∞ scale. A key advantage of the Bayesian approach is that it can be used to formally quantify the relative strength of the evidence for H₀ (“test the null”), in contrast to standard null-hypothesis significance tests (Wagenmakers et al., 2018; Bo et al., 2024). It also does not require the data analyst to choose what constitutes a trivial difference, unlike traditional equivalence tests (Hur et al., 2020b). BF₁₀ was interpreted using established benchmarks (van Doorn et al., 2021). Values <1 were interpreted as evidence of statistical equivalence (i.e., support for the null hypothesis), ranging from strong (BF₁₀ ≤ 0.10), to moderate (BF₁₀ = 0.10–0.33), to weak (BF₁₀ = 0.33–1). The reciprocal of BF₁₀ represents the relative likelihood of the null hypothesis (e.g., BF₁₀ = 0.10, H₀ is 10 times more likely than H₁).

At the behest of a reviewer, we also explored potential differences in BST and Ce function using activation estimates derived from the OSP model. Because the OSP model does not permit meaningful between-moment contrasts, we used paired Student's t tests to compare (1) sustained responses to uncertain-threat and (2) phasic responses to certain-threat anticipation, each relative to their second-level reference conditions (e.g., certain threat vs certain safety).

Brain-signature analyses

Standard fMRI analyses cannot address the momentary dynamics of threat-evoked emotions (Poldrack, 2011; Grogans et al., 2023). Here we used an independently trained and validated multivoxel pattern or “signature” of negative affect to covertly probe moment-by-moment fluctuations in anticipatory distress during the “countdown” period—something that would otherwise entail the imposition of a secondary rating task (e.g., using a continuous dial or randomized prompts), with unknown consequences for on-going emotional experience. Čeko, Wager, and colleagues used machine learning to develop a pattern of voxelwise weights predictive of the intensity of negative affect in unseen data (Čeko et al., 2022). They demonstrated that the signature is a sensitive indicator of distress elicited by a range of noxious experiences—including thermal and mechanical pain, unpleasant photographs, and aversive auditory stimuli—but unrelated to the intensity of feelings triggered by positive stimuli, indicating specificity. We computed the dot-product between the negative-affect signature and activation estimates derived for the present sample using the Convolved Blocks model, enabling us to generate signature responses (a probabilistic estimate of negative affect intensity) for every combination of threat certainty (certain, uncertain), block (1, 2, 3), and participant. We used one-sample Student's t tests to confirm that the signature, which was trained using activation estimates time locked to the presentation of aversive stimuli, is nominally sensitive to the anticipation of threat encounters (p < 0.05, uncorrected). For hypothesis testing, we used a standard 2 (Threat Certainty: Certain, Uncertain) × 3 (Convolved Block: 1, 2, 3) repeated-measures GLM (Huynh–Feldt correction) to assess dynamic fluctuations in signature-estimated distress across threat contexts. Significant GLM effects were decomposed using polynomial contrasts, as detailed in the prior section.