Bed nucleus of the stria terminalis connectivity during food cue and taste processing under stress

Domain: neuroscience
Published: 2026-04-08
Editorial outcome: major_revision

Abstract
Across species, stress drives alterations in feeding behaviour, including heightened food-seeking and the overconsumption of palatable foods. The bed nucleus of the stria terminalis (BNST) acts as a neural hub linking stress and reward circuits. However, its role in human food cue and taste processing under stress remains unclear. Here, using 7-Tesla fMRI, dynamic causal modelling and in-scanner taste delivery, we demonstrate that beverage cues and taste receipt under stress modulate BNST effective connectivity. Forty-eight participants were presented a palatable and neutral beverage cue before receiving the corresponding beverages under low- and high-stress conditions. Beverage cues under stress downregulate BNST effective connectivity to the nucleus accumbens, orbitofrontal cortex, and dorsal mid-insula (dmINS), with the strength of cue-related downregulation of BNST-to-orbitofrontal effective connectivity related to changes in participant’s subjective stress. In addition, taste receipt under stress downregulates effective connectivity from the dmINS to BNST. These findings provide evidence that stress-related food cue and taste processing leads to modulation of BNST reward circuitry, providing a mechanistic account for how stress modulates responsivity to food cues.

Introduction
Many individuals report increasing their food intake in response to stress1,2, and stress is associated with eating pathology, particularly binge eating3. Though the neural mechanisms remain unclear, stress enhances responsivity to food cues, motivates food-seeking behaviour4,5, and shifts preferences towards palatable food6 irrespective of weight status7. Furthermore, stress can increase cue-triggered craving of sweet rewards in humans8 and lower taste intensity9,10, which may in turn increases sweet sensitivity, an individual difference associated with greater liking of sweet foods11. Rewarding high-calorie food cues, which are abundant in the present food environment, can also trigger appetitive drive in the absence of satiety12, and promote the overconsumption of food under stress13. Thus, in this way, the stressed brain appears to be biased towards an increased desire to eat, coupled with a reduced ability to inhibit feeding behaviour, which is precipitated by the presence of rewarding food cues1.

The bed nucleus of the stria terminalis (BNST) serve as a primary hub for integrating internal bodily cues that guide motivated behaviour—such as feeding—in response to threats and rewards14,15. The BNST is anatomically connected to the amygdala, lateral hypothalamus, striatum, insula and orbitofrontal cortex (OFC)16–18, regions implicated in stress, eating and reward responses, respectively. In sated mice, optogenetic activation of the mostly GABAergic (inhibitory) neurons projecting from the BNST to the lateral hypothalamus induces rapid consumption of palatable food19,20. Human BNST connectivity under food cue and taste processing during stress, however, has not been well characterised.

To date, human BNST research has focused primarily on its role in stress, anxiety and drug reward, including stress-induced drug-seeking behaviour21. In active smokers, enhanced cue responsiveness to smoking videos is predicted by reduced BNST, nucleus accumbens (NAcc) and insula activation during acute psychosocial stress22. Other activation-based studies implicate the NAcc, insula and OFC in heightened food cue responsiveness under stress23,24, and in the encoding of food taste25–27. Evidence from diffusion-weighted MRI, resting-state fMRI and post-mortem studies have confirmed bidirectional white matter tracts between the BNST, NAcc, insula and OFC in humans28–32, such that these regions form an anatomical network. Within this network, Banasikowski and Hawken14 have suggested the BNST may play a role in ‘valence surveillance’33, monitoring and integrating information about the emotional value of stimuli by relaying hunger and satiety signals from the hypothalamus to the OFC, which is involved in goal-directed decision-making.

Mapping the BNST’s function in humans at standard MRI field strengths is challenging due to its small size35. By virtue of its superior spatial resolution, ultra-high-field 7-Tesla (7T) fMRI overcomes this barrier, permitting the resolution of neural activity in deep brain structures, including the BNST32. Using such approaches, progress has been made with regard to characterising the food cue-related activity of large-scale brain networks36–41. However, understanding food cue-related interactions within brain networks has been hampered by the field’s reliance on correlational measures (functional connectivity, for instance), rather than methods that can infer effective connectivity—the signed, directed influence that one brain region exerts over another42, where a positive (negative) connection reflects a net excitatory (inhibitory) influence at the level of neuronal populations. Although not necessarily superior to simpler functional connectivity approaches in all contexts, dynamic causal modelling (DCM) offers a biophysically grounded framework for inferring effective connectivity, and has proven valuable in addressing such questions43.

Here, we leveraged 7T fMRI, DCM and controlled taste delivery via a gustometer to map the effective connectivity between the BNST and key regions implicated in food cue and taste processing: the NAcc, mid-insula (dmINS), and OFC. We hypothesised that food cue and taste processing under stress would downregulate (decrease the values of) BNST connectivity with other network nodes. We further hypothesised that effective connectivity within the BNST network would be influenced by the reward value of beverage cues and taste. Our findings show that beverage cues under stress have a down-regulatory effect on BNST connectivity to the NAcc, OFC and dmINS. In addition, the strength of the high-stress cue-related downregulation of BNST-to-OFC effective connectivity could predict individual-participant changes in subjective stress, providing evidence of BNST-connectivity–behaviour relationships in the context of food cue processing.

Results
Beverage task and stress induction
All participants (n = 48) completed two runs of a block-design beverage task: one following low-stress induction and another following high-stress induction (Fig. 1; detailed in Methods). During each task run, participants were shown images of either a glass of chocolate milk or water. Following each cue presentation, participants received the corresponding beverage inside the scanner via a gustometer.

Self-reported stress levels were higher after the high-stress induction (mean 5.04 ± 1.58 [s.d.]) compared to the low-stress induction (mean 2.27 ± 1.48, t = 12.41, p = 1.94 × 10-16, two-sided, n = 48; Fig. 2a). There was a main effect of stress (F(1, 94) = 15.01, p = .0002, two-sided, n = 48) and beverage (F(1,94) = 38.91, p = 1.26 × 10-8, two-sided, n = 48) on mean pleasantness ratings: chocolate milk was rated more pleasant than water after both high-stress (mean 5.26 ± 1.23 versus mean 3.88 ±1.26, p = 2.85 × 10-8) and low-stress induction tasks (mean 5.61 ± 1.12 versus mean 4.29 ± 1.04, p = 9.00 × 10-8), two-sided for both; Fig. 2b, Supporting Information). Analysis of cue-level pleasantness ratings revealed a main effect of cue exposure (F(5, 940) = 10.17, p = 1.66 × 10-9, two-sided), an interaction between cue exposure and stress (F(5, 940) = 2.64, p = .02, two-sided), and no significant interaction between stress, beverage type and cue exposure (F(5, 940) = 1.007, p = 0.41, two-sided; Fig. 2c).

Bonferroni-corrected post hoc pairwise comparisons using paired t-tests showed that chocolate milk was rated more pleasant than water at all cue exposures during high stress (cue 1: p = .0004; cue 2: p = .0019; cue 3: p = 1.85 × 10-5; cue 4: p = 1.21 × 10-5; cue 5: p = .0002; cue 6: p = 6.33 × 10-5) and at all cue exposures except cue 1 during low stress (cue 2: p = .0002; cue 3: p = 1.85 × 10-5; cue 4: p = .0006; cue 5: p = 4.22 × 10-5; cue 6: p = 1.85 × 10-5), two-sided, n = 48 for all. These results suggest that (1) the stress manipulation was successful in increasing participants’ self-reported stress levels, and (2) participants rated the rewarding beverage (chocolate milk) as more pleasant than the neutral beverage (water) throughout the beverage task.

Physiological stress was primarily indexed using tonic electrodermal activity (EDA; detailed in Methods), complemented by post hoc analyses of heart rate variability (HRV). EDA values were not significantly different after high-stress induction compared to the low-stress induction (t(27) = -0.312, p = .758, two-sided, n = 28, mean difference = -0.006 ± 0.11; see Fig. S1). Heart rate variability, quantified as the root mean square of successive differences (RMSSD), was lower at the beginning of the high-stress run than at the beginning of the low-stress run (t(43) = -2.25, p = 0.029, two-sided), consistent with greater physiological stress during the high-stress condition. We also used sliding windows (10-second windows stepped every 2 seconds) to analyse the complete time course of RMSSD across both task runs. Data were segmented into early, middle, and late task periods (first, middle, and final 20 seconds) to evaluate trends across time. A within-participants repeated-measures ANOVA found a significant condition × time interaction (F(2, 86) = 3.20, p = 0.045, two-sided), indicating that RMSSD evolved differently over time between the low- and high-stress conditions (Fig. S2). Uncorrected post hoc comparisons showed that RMSSD increased significantly more from early-to-middle in the low-stress condition compared to the high-stress condition (t(43) = 2.25, p = 0.029, two-sided). Likewise, RMSSD increased in the low-stress condition but decreased in the high-stress condition when comparing early-to-late windows (t(43) = 2.03, p = 0.048, two-sided). Additional short-term HRV indices (pNN20 and pNN50; the proportion of successive inter-beat intervals differing by >20 ms and >50 ms, respectively) and a non-linear Poincaré-derived measure (SD1/SD2; the ratio of short- to long-term variability) showed a consistent pattern indicative of greater physiological arousal during the high-stress run. Mean heart rate averaged across the entire run did not differ between conditions (t(47) = −0.37, p = 0.710, two-sided). Full results are reported in the Supplementary Information.

Brain activity during beverage cue presentation and taste receipt
A whole-brain, group-level general linear model (GLM) revealed significant activation in response to beverage cue presentation (‘CUE’), with robust effects present in visuospatial and attentional regions including the lingual gyrus, angular gyrus and middle frontal gyrus. Additional significant clusters were found in the dorsolateral prefrontal cortex, mid-cingulate cortex, middle temporal gyrus, and middle frontal gyrus, with the OFC and dmINS showing activation, overlapping with temporal pole clusters. Subcortical activations included the putamen, thalamus, hippocampus, ventral striatum (NAcc) and amygdala extending to the BNST, with ancillary recruitment of the substantia nigra and ventral tegmental area (VTA; pFDR < .05, one-sided, n = 48, voxel-wise false-discovery-rate [FDR] corrected; Fig. 2d).

A second group-level GLM analysis revealed significant activation in response to beverage taste receipt (‘TASTE’), with the strongest effects in sensorimotor and gustatory regions, including the postcentral and precentral gyri, supplementary motor area, cerebellum, rolandic operculum and dmINS. Subcortical engagement was notable in the putamen and thalamus. The OFC showed modest engagement, alongside discernible activation in the amygdala extending to the BNST, and ancillary recruitment of the NAcc (pFDR < .05, one-sided, n = 48, voxel-wise FDR corrected; Fig. 2e, for complete results, Tables S1–S2). Shared activation across cue and taste conditions for the DCM volumes of interest (VOIs) are depicted in Fig. S3.

Bed nucleus of the stria terminalis effective connectivity modulation under stress
As suggested by published guidelines for DCM44,45, the selection of VOIs for analysis was informed by GLM activation results to ensure that VOIs contained task-relevant activity. The dynamic causal model network structure comprised four regions: the BNST, NAcc, dmINS and OFC (Table 1). A priori, we permitted interactions to conform to a star-shaped connectivity architecture, with the BNST as the central hub (Fig. 3). To examine how beverage cue presentation and taste receipt modulates BNST connectivity under stress, for each participant, four separate dynamic causal models were specified and inverted using concatenated time-series data from the low- and high-stress runs. Each model assessed the distinct modulatory effect of: (1) high-stress cues: the modulatory effect of all beverage cue presentations under high stress, (2) high-stress chocolate milk cues: the modulatory effect of rewarding (chocolate milk) beverage cue presentation under high stress, (3) high-stress taste: the modulatory effect of all beverage taste receipt under high stress, (4) high-stress chocolate milk taste: the modulatory effect of rewarding (chocolate milk) beverage taste receipt under high stress. For the ‘CUE’ models (Models 1 and 2), the driving input (to all regions) was specified as all beverage cue onsets across both runs. For the ‘TASTE’ models (Models 3 and 4), the driving input was specified as all beverage taste receipt onsets across both runs. The modulatory inputs on the strength of each connection to and from the BNST were condition-specific: (1) all beverage cue onsets under high-stress, (2) only chocolate milk cue onsets under high-stress, (3) all beverage taste receipt onsets under high-stress, (4) only chocolate milk taste receipt onsets under high-stress. Four participants were excluded from the ‘CUE’ DCM analyses (n = 44) and one from the ‘TASTE’ DCM analyses (n = 47) as they did not survive the BNST VOI time-series extraction (for details, see Methods). Results from additional models testing the effects of the low-stress and water conditions are included in the Supporting Information (Tables S3–S4). Results from additional over-complete models—specified identically to the main analyses except that the modulatory inputs were coded to represent all four conditions simultaneously—are also reported in the Supporting Information (Table S5, Fig. S4). These analyses suggest that modulation is primarily stress-driven, with minimal evidence for stimulus-specific (water versus chocolate milk) effects.

Parametric empirical Bayes (PEB) was used to summarise group-level effects from DCM and to assess the association between the cross-condition change in self-reported stress ratings—change in subjective stress—and the modulation of effective connectivity. We found evidence to suggest that high-stress beverage cue presentation downregulated efferent connectivity from the BNST to the NAcc, OFC, and dmINS, as well as efferent connectivity from the dmINS to the BNST (Fig. 4a). When focusing on chocolate milk cues (Fig. 4b), we found that cue presentation under high stress downregulated efferent connectivity from the BNST to dmINS, and upregulated (increased the value of) efferent connectivity from the NAcc to BNST. By contrast, focusing on taste receipt under stress, we found that taste receipt downregulated dmINS-to-BNST effective connectivity (Fig. 4c); however, this effect did not exceed our imposed posterior probability threshold when considering chocolate milk and water taste receipt separately.

When focusing on both high-stress beverage cues and treating change in subjective stress as a covariate, we found that downregulation of BNST-to-OFC effective connectivity was associated with change in subjective stress such that individuals with greater change in subjective stress, showed less downregulation (Fig. 4d). Using leave-one-out cross-validation, we showed that this effect was robust: across folds, the group level model of BNST-to-OFC effective connectivity could predict the change in subjective stress of left-out participants (r = .28; p = 0.031, two-sided; Fig. 4e–f). Posterior probabilities and modulatory connections shown in Fig. 4a–d, are reported in Table 2.

Methods
Participants
This study was approved by the Human Research Ethics Committee of The University of Melbourne (HREC ID 22347) and conducted in accordance with all relevant ethical regulations. Fifty-seven healthy adults were recruited into the study via online advertisements. Participants were eligible if they (1) were 18 to 40 years old, (2) had no history of eating disorder symptomology and/or current diagnosed mental health disorder, (3) were fluent in verbal and written English, (4) had no major hearing impairment or damage, (5) were able to consume dairy products, (6) had no contraindications to MRI (for example, pregnancy, claustrophobia or metal present in their body they cannot remove), and (7) had no medical or neurological condition for which they were on medication. Participants were screened for mental health disorders using the DSM-IV/ICD-10 version of the Mini-International Neuropsychiatric Interview (MINI)82. Forty-eight participants (mean age 26.2 ± 5.3 years; 24 male, 24 female; self-reported) successfully completed the scanning protocol and were included in GLM analyses (further demographics, Table S6). Forty-four participants were included for the ‘CUE’ DCM analyses and 47 for the ‘TASTE’ DCM analyses. Participants gave their written informed consent before commencing the study at the Melbourne Brain Centre Imaging Unit (The University of Melbourne, Parkville).

Beverage task and stress induction
To investigate the effect of food cue and taste processing under stress on neural connectivity, all participants completed two runs of a block-design beverage task that utilised an MRI-compatible gustometer (Multistimulator OG001, Burghart Messtechnik, Pinneberg, DE) and PsychoPy83 v2024.2.4. Prior to scanning, all participants received training on how to complete the task, were shown the beverage cues, and explicitly told they would receive the matching beverages in the MRI.

Prior to each beverage run, participants completed a structured stress induction task. The first induction (low-stress condition) involved an easy, self-paced serial subtraction task with no evaluative feedback, while the second (high-stress condition) incorporated a difficult, time-pressured serial subtraction task with negative social evaluative feedback, modelled on the Montreal Imaging Stress Task (MIST)84–87. In both conditions, participants viewed an on-screen digital clock (a simple rotating arm which remained white in colour for the first minute and turned red for the second minute). There is a strong precedent for using brief, in-scanner math-related tasks to induce stress in the context of fMRI84–87. Further information on the stress induction procedure is available in the SI Methods.

During each 7-minute beverage task run, participants were shown an image (‘CUE’ condition) of a glass containing either chocolate milk (rewarding cue) or water (neutral cue) on an MRI-compatible screen for 5 s (Fig. 1). Following each image, participants received 1 ml of the corresponding beverage (‘TASTE’ condition), whilst the word “TASTE” was displayed for 5 s to prompt participants to taste and swallow the beverage. Chocolate milk (Table S7) or a tasteless water solution (25 mM potassium chloride, 2.5 mM sodium carbonate dissolved in water) were delivered via a 3D-printed mouthpiece88 (SI Methods). Following beverage receipt, participants had 5 s to rate—using an MRI-compatible button box—how pleasant the beverage tasted on a 7-point Likert scale, with 1 being very unpleasant and 7 being very pleasant. After each rating, “RINSE” appeared on the screen for 5 s, at which point participants received 1.5 ml of the tasteless solution via the mouthpiece to use as a rinse. Each run consisted of 12 blocks (6 blocks of chocolate milk, 6 blocks of tasteless water solution). The order of chocolate milk and water beverage presentation and corresponding delivery was pseudorandomised. Prior to the start of each run, participants rated their hunger and stress levels on a 7-point Likert scale with 1 being not hungry (stressed) and 7 being very hungry (stressed) (Fig. 1). All participants provided responses for beverage, hunger and stress ratings across both beverage task runs.

Behavioural analyses
Self-reported hunger, stress and chocolate milk pleasantness ratings were tested for normality using the Shapiro-Wilk method. Hunger and stress levels after low- versus high-stress induction tasks were analysed using two-sided paired t-tests. To assess overall beverage pleasantness, mean ratings of chocolate milk versus water were analysed with a two-way repeated measures ANOVA (stress × beverage), followed by Bonferroni-corrected paired t-tests for post hoc comparison. To examine cue-specific pleasantness dynamics, ratings of each beverage during each cue (taste) block of the low- and high-stress runs were analysed with a three-way repeated measures ANOVA (beverage × stress × cue) with Bonferroni-corrected paired t-tests at each cue exposure. All post-hoc tests were two-sided and n = 48 for all comparisons. Data were analysed and visualised in GraphPad PRISM v10.2.2 (GraphPad Software, Boston, MA, US). Demographic data were analysed using SPSS v29.0.2.0 (IBM Corp., Armonk, NY, USA, 2023).

Physiological measures
EDA data were recorded during both low- and high-stress induction tasks. These data were recorded using MRI-compatible finger electrodes (Ag/AgCl) coated in conductance gel (0.5% saline) and placed on the index and middle finger intermediate phalanges of participants’ left hands. The signal was amplified and sampled at 1000 Hz using PowerLab version 8.0 (AD Instruments).

EDA data underwent manual review for quality by two independent raters (E.G.-H. and M.D.G.) both before and after preprocessing, which included filtering and downsampling. Filtering was performed using a first-order Butterworth bandpass filter with cutoff frequencies set at 0.0159 Hz and 5 Hz. Following filtering, the signals were downsampled to a target rate of 10 Hz. Raters classified files as acceptable, discardable, or requiring the imputation of a missing event marker. Files exhibiting pervasive artifacts, missing both start and end markers, or flat responses (indicative of non-responses) were excluded from further analysis. The inter-rater reliability for participant inclusion decisions was 78.3%, and after resolving discrepancies between raters, data from 28 participants with recordings for both stress-induction tasks remained.

To infer tonic sympathetic arousal from the spontaneous skin conductance fluctuations occurring during both stress-induction tasks, a matching pursuit algorithm89 was applied using the psychophysiological modelling (PsPM) toolbox90,91 (version 6.1.2). The algorithm was configured to model sympathetic arousal without applying additional filtering beyond the initial preprocessing, and a threshold of 0.1 µS was set for peak detection. A paired two-sided t-test was used to compare sympathetic arousal between runs.

The RMSSD in pulse signals extracted using the PhysIO toolbox were used as a measure of HRV during the low- and high-stress conditions. Pulse oximetry was only collected and synchronised during functional image acquisition, meaning HRV estimates could only be derived from the beverage task runs. In accordance with previous research92,93, a window of 20 seconds with 2-second steps was used to assess RMSSD differences at the beginning of each run. As stress responses unfold dynamically within a run, collapsing RMSSD across the entire task can obscure periods of heightened arousal and subsequent recovery, whereas computing RMSSD in shorter sliding windows provides the sensitivity needed to detect these transient physiological shifts. Indeed, mean RMSSD did not differ when examining the entirety of low- and high-stress runs (t(43) = −0.60, p = 0.55, two-sided), whereas significant differences did emerge using a sliding-windows approach.

RMSSD values were extracted continuously across the task period and subsequently segmented into early, middle, and late phases, defined as the first, middle, and final 20 seconds of the task, respectively. These values were entered into a within-participants repeated-measures ANOVA with condition (low-stress, high-stress) and time (early, middle, late) as factors. This design follows recommendations to model HRV as a time-varying signal, allowing condition effects to be tested while capturing the distinct phases of autonomic reactivity within each run94,95. Post hoc comparisons between time points within each condition were conducted using paired t-tests to assess directional trends over the course of the task. In addition, successive inter-beat intervals were converted into difference vectors, from which the proportion of successive differences exceeding 20 ms (pNN20) and 50 ms (pNN50) was calculated. The nonlinear Poincaré index (SD1/SD2) was derived from the standard deviations of instantaneous beat-to-beat changes (SD1) and longer-term variability (SD2). Mean values from low- and high-stress runs were compared using paired t-tests.

Functional MRI acquisition and preprocessing
MRI acquisition was conducted at the Melbourne Brain Centre Imaging Unit, The University of Melbourne, using a 7T research scanner (Siemens Healthcare, Erlangen, DE) fitted with a 32Rx/1Tx-channel head-coil (Nova Medical Inc., Wilmington, MA, US). The functional sequence (beverage task) used a gradient-echo echo-planar imaging (GE-EPI) acquisition with multiband acceleration (factor 6) and GRAPPA (factor 2) in the steady state (repetition time 800 ms; echo time 22.2 ms; flip angle 45°; field of view 20.8 cm; slice thickness 1.6 mm, no gap; matrix 130 × 130; 84 interleaved axial slices aligned to the anterior–posterior commissure line)96. The run duration was 7 min, yielding 525 whole-brain EPI volumes per session. A high-resolution T1-weighted, multi-echo MP2RAGE dataset was acquired at the start of the MRI protocol (repetition time 4500 ms; echo times 2.21/4.21/6.15/8.14 ms; inversion times 700/2700 ms; flip angles 6/7°; field of view 24 cm; slice thickness 0.75 mm; matrix 320 × 320; 224 axial interleaved slices)97. The denoised, combined-echo 3D MP2RAGE image was used to facilitate functional image co-registration98. Participants’ heads were stabilised using foam-padded inserts. Cardiac and respiratory signals were sampled at 50 Hz using a Siemens pulse oximeter and respiratory belt and were used for physiological noise correction.

Structural and functional data were pre-processed using the statistical parametric mapping toolbox (SPM12; v7771, Wellcome Trust Centre for Neuroimaging, London, UK) in a MATLAB 2023a environment (The MathWorks Inc., Natick, MA, US). To correct for motion artefacts, individual participants’ time-series were realigned to the mean image and individualised motion regressors created using the Motion Fingerprint Toolbox99. Anatomical images were co-registered to their mean functional image, segmented and normalised according to the International Consortium of Brain Mapping template using the unified segmentation plus diffeomorphic anatomical registration through exponentiated lie algebra (DARTEL)100. Images were then smoothed with a 3.2 mm full-width-at-half-maximum Gaussian kernel to preserve spatial specificity.

Physiological noise was included as a regressor by using the PhysIO toolbox at the first level101. This approach uses physiological recordings to apply noise correction to fMRI sequences and has previously been shown to benefit the temporal signal-to-noise-ratio of blood-oxygen-level-dependent (BOLD)102,103. Specifically, we used retrospective image-based correction104, to model periodic cardiac and respiratory contributions to the BOLD signal using canonical respiratory105 and cardiac response functions106. Individualised DARTEL-derived tissue masks were used for anatomical component correction (aCompCor107) by extracting nuisance components from white matter and cerebrospinal fluid signals.

General linear models of beverage cue and taste processing
First-level (single-participant) GLMs were specified separately for each beverage-task run using the preprocessed fMRI time series and participant-specific nuisance regressors (head motion and physiological noise). Regressors of interest comprised the onsets of beverage cue presentation (‘CUE’) and taste receipt (‘TASTE’) for chocolate milk and water, each convolved with the canonical haemodynamic response function. All remaining task events (for example, pleasantness ratings and the rinse period) were modelled as regressors of no interest. A 128-s high-pass filter was applied to remove low-frequency drift. Serial autocorrelation was modelled using SPM’s FAST approach108, and models were estimated in SPM using prewhitening and generalised least squares.

For group-level inference and to identify candidate regions for subsequent DCM analyses, subject-wise first-level contrast images were entered into a second-level GLM implemented as a voxelwise paired t-test, testing whether the mean within-participant contrast differed from zero across participants. Separate group-level analyses were performed for ‘CUE’ (all cue events collapsed) and ‘TASTE’ (all taste events collapsed). Whole-brain results were thresholded using voxel-wise FDR correction (pFDR < 0.05) with a minimum cluster extent of 10 voxels.

Dynamic causal modelling
DCM is used to infer directed influences between neuronal populations from observed BOLD time series data. A model is inverted (fitted using a variational Laplace, a Bayesian optimisation technique) to obtain the model structure and connectivity parameters that maximise parsimony: balancing accuracy and complexity43. This modelling framework can be utilised to infer the extent to which a manipulated input—an experimental task, for example—upregulates (positively modulates) or downregulates (negatively modulates) the intrinsic excitatory (positive) or inhibitory (negative) source-to-target communication between neuronal populations (nodes). Such continuous-time modulation can be framed in terms of rates (measured in Hz), whereby a task modulates the degree to which activity of an efferent node affects the continuous rate at which the state of afferent node changes. Previous evidence has demonstrated that 7T fMRI furnishes more reliable effective connectivity inferences compared to lower MRI field strengths109.

Region selection and time-series extraction
Four regions showing significant activation during ‘CUE’ and ‘TASTE’ were included in dynamic causal models: the BNST, the NAcc, the OFC and the dmINS. The VOIs for each region were extracted at the subject level as described in published guidelines44. For the dmINS and OFC, the first eigenvariate across all activated voxels within a 4 mm radial sphere around the subject-specific maxima (p < .05, uncorrected), but no more than 8 mm from the group maxima, were extracted. Time-series extraction for the BNST was run using a consensus mask, created by overlaying the manually segmented BNST mask from Theiss and colleagues110, onto Pauli and colleagues’ extended amygdala mask111. Both masks were binarised and combined to create a new mask containing only voxels common to the original masks. The NAcc mask was generated using the 1-mm3 resolution atlas from automated anatomical labelling 3 (AAL3112) toolbox. The OFC and dmINS were modelled unilaterally, using the strongest task-evoked peaks from the GLM, in line with DCM guidelines suggesting prioritisation of the dominant hemispheric activation (which in this sample, were located in the right hemisphere). The BNST and NAcc were modelled bilaterally given their small size, proximity to the midline, and known bilateral function. Four participants were excluded from the ‘CUE’ DCM analyses (n = 44) and one from the ‘TASTE’ DCM analyses (n = 47) during regional time-series extraction as no significant BNST voxels passed the significance threshold. Individual-level VOI peak locations are depicted in Fig. S5.

To model stress-related differences within a single dynamic causal model per participant, time series from the low- and high-stress runs were concatenated using spm_fmri_concatenate. Regional time series were extracted as GLM-adjusted eigenvariates: nuisance effects were regressed out during VOI extraction by adjusting the time series for regressors not included in the effects-of-interest F-contrast. The resulting residualised regional signals (and associated observation noise) were carried forward into DCM.

Model specification and inversion
Based on previous anatomical literature16–18, DCM structures were specified such that they assumed bidirectional intrinsic connections between the BNST and the NAcc, OFC and dmINS (Fig. 3). At this first level, models were inverted using concatenated time series data from both low- and high-stress task runs (under standard, default priors). The cue models included all high-stress and low-stress ‘CUE’ blocks (n = 44) as the first input condition, high-stress ‘CUE’ blocks as the second input condition (high-stress cue) and high-stress rewarding ‘CUE’ blocks as the third input condition (high-stress chocolate milk cues). This structure was repeated for the taste models and included all high-stress and low-stress ‘TASTE’ blocks (n = 47). The modulatory effect of food cues, taste, rewarding cues (chocolate milk cues) and rewarding taste (chocolate milk taste) under high stress were tested on all pathways to and from the BNST.

Each model assessed a distinct modulatory effect on BNST effective connectivity: (1) high-stress cues: the modulatory effect of all beverage cue presentations under high stress, (2) high-stress chocolate milk cues: the modulatory effect of rewarding beverage cue presentation under high stress, (3) high-stress taste: the modulatory effect of all beverage taste receipt under high stress, (4) high-stress chocolate milk taste: the modulatory effect of rewarding beverage taste receipt under high stress. For the cue models (Models 1 and 2), the driving input was specified as all beverage cue onsets across both runs. For the taste models (Models 3 and 4), the driving input was specified as all beverage taste receipt onsets across both runs.

Using four separate dynamic causal models instead of a single model with multiple modulatory inputs preserves the theoretical and experimental distinctions between stress conditions and cue- and taste-related processes, preventing interactions that could obscure meaningful differences43. It also reduces model complexity (overparameterisation), and helps to ensure that connectivity changes are clearly attributable to specific conditions43.

Parametric empirical Bayes
Inversion of dynamic causal models yielded posterior distributions over effective connectivity parameters. These posterior distributions for the modulatory effects served as data in PEB models wherein subject-level effective connectivity represents random deviations from an average (group-level) effective connectivity network. One PEB model (Fig. 4d) included the effect of change in subjective stress, quantified as the change (delta) from low- to high-stress induction45 (specifically, the self-reported stress score after high-stress induction minus the self-reported stress score after low-stress induction). After inverting each PEB model, Bayesian model reduction and a greedy search were utilised to evaluate reduced models by omitting connections. Plausible reduced models were averaged using a Bayesian model averaging scheme (per the spm_dcm_peb_bmc routine). Evidence of a nonzero group effect was inferred when the two-sided posterior probability exceeded 0.75 (positive evidence113).

Leave-one-out cross-validation
Leave-one-out cross-validation was used to test the predictive accuracy of group effect sizes on participants’ individual changes in subjective stress, using the spm_dcm_loo routine. Focusing on the high-stress cue-related modulation of the BNST-to-OFC connection, a PEB model was iteratively inverted for all but one participant to test whether the modulation of effective connectivity parameters was robust enough to predict the change in subjective stress of the left-out participant. The reliability of this prediction was indicated by the correlation between predicted and observed values. An additional subject-level analysis (Fig. S6) showed that individual differences in stress-evoked BNST-to-OFC modulation (encoded in dynamic causal models) was positively associated with changes in subjective stress, providing a complementary, non-PEB perspective on the same effect.

In silico validation of dynamic causal models
The accuracy of DCM can be affected by factors such as signal quality. To address this issue, we conducted an in silico analysis of the main cue- and taste-related dynamic causal models to assess robustness of parameter recovery to estimated noise. These analyses involved three steps: (1) estimation of subject- and region-specific signal-to-noise ratios (SNRs) from empirical BOLD time-series data; (2) simulation of synthetic BOLD time-series using inferred model parameters—including hemodynamic delay parameters—that was combined with region-specific noise reflecting the estimated SNRs; and (3) model inversion using both DCM and PEB procedures on the simulated data. We then evaluated parameter recovery and the correspondence between second-level estimates derived from empirical versus simulated data (SI Results; Table S8; Figs. S7–S8).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript by Guerrero-Hreins et al. seeks to address the highly relevant question of whether stress modulates BNST connectivity during taste processing. The study has several strengths, such as using a high temporal resolution 7T sequence that provides a strong resource for dynamic causal modelling, an innovative approach based on the experimentally demanding use of an MR-compatible gustometer, and a sufficiently large sample. Despite these strengths, there are multiple major issues that require substantial revision before the manuscript can provide the aspired answer. Illustratively, I would expect far more dedication to the right use of anatomical regions if this is one of the key contributions of the paper, and I have doubts whether the current definition is justified.\nMajor issues\n1. Although the introduction is written in an accessible manner, the statement concerning the rationale for DCM should be rephrased. I understand that the authors want to argue that this method is better than (undirected) FC, but it is only superior to FC if certain assumptions are met (and ideally evaluated). For example, it is known that a differential hemodynamic delay can lead to spurious inferences about the directionality of the effects. Likewise, differences in tSNR or the seeds can contribute to incorrect inferences. Without additional control analyses, it is difficult to take the (strongly worded) inferences for granted.\n2. The phrasing of the (clearly) non-significant EDA effects is odd. If the numbers are correct, it is not even close to a directional effect. In my opinion, this should be clearly stated. Since EDA failed to capture stress-related physiological effects, the authors may consider using alternatives (e.g., pulse rate as in (Kühnel et al., 2022; Erhart et al., 2024))\n3. In general, the use of causal language for functions of brain regions (“OFC is responsible for”) is not justified unless there is very clear evidence from lesion studies. I strongly suggest rephrasing these sections to tone down these statements.\n4. The labels for the seed regions are highly misleading in two cases. First, the selected seed is not in the anterior insula, but close to the dorsal mid insula, likely including the primary ingestive cortex. This is plausible given the design of the task. Still, it weakens the argument that direct structural connections to the BNST exist (they are mostly to the anterior part of the insula, which is also involved in taste processing, but likely more concerning expectations). Second, the BNST highlighted in the figures does not correspond well to the BNST. The reference to the Pauli atlas is also misleading in that sense because they only provide a mask for the extended amygdala. This includes the BNST, but also other parts, and it seems as if the authors have used the entire mask as search volume, leading to a situation where the group-level peaks are located well outside the BNST. Also, the provided coordinate numbers do not correspond to the figure, so I guess a minus sign is missing for the z coordinate (otherwise, they would be placed even further away from the targets). In my opinion, this is a major problem since the paper is focused on the function of the BNST while the current analysis applies more to the extended amygdala, which is much more difficult to connect to the preclinical work.\n5. It is not clear to me what the evidence is that BNST input would inhibit goal-directed behavior as there was no association with behavior that would support this (should be rephrased). In general, the authors do not make much use of the behavioral ratings to strengthen their behavioral claims in the discussion. I think they should provide additional (exploratory) analyses with changes in stress and liking or edit the (far too strongly worded) language in the discussion.\n6. The discussion should be rebalanced to better reflect which findings speak against the (posited) specificity. For example, isn’t it surprising in light of the hypothesis that beverage taste receipt under stress elicited similar inhibitory effects and was not specific to the milkshake? How would the authors rule out that they are reporting a less functionally distinct process here? Water is essentially also a control condition for everything that is simply driven by the timing of the task elements, and I do not think this interpretation was given enough consideration in the report. This concern is even more important given the lack of randomization of the block order.\n7. Please restructure the information so that the paper can be read linearly. For example, it is only later revealed why the authors excluded 3 participants for the analysis of the cue responses? Were the exclusion criteria prespecified, and would the TASTE DCM with the same 45 participants give the same results? Wouldn’t this approach inflate the model estimates if people without a clear activation in VOIs are selectively excluded (was it due to the BNST or various regions)? Couldn’t the subject-specific peak from the TASTE analysis be used instead?\n8. Perhaps I missed it, but I do not see any attempt to account for the potential time-related effect of the preceding stressor. I understand that it is difficult to dissociate linear trends from signal drift in the MR signal, but the beverage runs are quite long compared to the preceding stress exposure. Was the available physiological data (e.g., changes in pulse rate) explored to evaluate how long the after-effects of stress last?\n9. The sequence is optimized for acquisition speed, and this may improve the effective connectivity estimates. However, in my experience, the use of GRAPPA and a multiband factor >4 can substantially reduce the SNR in deeper brain regions. Did the authors calculate and evaluate tSNR in their target regions, since differences in tSNR (e.g., between the insula and the BNST) could also affect the model estimates of effective connectivity (since BNST signal will likely fluctuate much more).\nReferences\nErhart M, Kühnel A, Brückl TM, Fietz J, Pöhlchen D, Yusupov N, Elbau I, Czisch M, Binder EB, Sämann PG, Spoormaker VI (2024) Capturing inter-individual variability in stress dynamics with heart rate traces reveals activity in the bilateral hippocampus, amygdala, and insula. Imaging Neuroscience 2:1-19.\nKühnel A, Czisch M, Sämann PG, BeCOME working group, Binder EB, Kroemer NB (2022) Spatio-temporal dynamics of stress-induced network reconfigurations reflect negative affectivity. Biol Psychiatry 92:158-169.", "reviewer_2": "The authors present a fascinating paper related to the role of the bed nucleus of the stria terminalis (BNST) in stress-related feeding. The authors used a form of psychosocial stress to induce stress in participants, along with CUE and TASTE conditions related to the receipt of water or chocolate milk, respectively. The authors also used dynamic causal modeling to better determine direction of inhibition and excitation within a network comprised of the BNST, nucleus accumbens (NAcc), orbitofrontal cortex (OFC), and anterior insula (aInsula). The strengths of the paper are rationale based on good evidence from rodent models, strong experimental task, and generally excellent methods including covarying for physiological signal and leave one out analyses. While this is a great paper, there are several major and minor concerns related to the methodology and discussion of the data that reduce enthusiasm.\nMajor concerns:\n- The major concern is the definition of the brain regions used in the analyses. The area of activation listed as BNST in Figure 1 appears to be dorsal amygdala and not BNST. For the coordinates shown in Table 1, maybe there is an error for the BNST peak (should be z = -11?); the peak coordinates for the NAcc appears to be in the very most anterior portion, very close to white matter; the peak coordinates for the anterior insula appear to be in the dorsal middle or posterior insula (see 10.1093/scan/nss066). It would be helpful to overlay the areas used in the DCM analysis on Figure 1. For the enlarged parts of Figure 1, please show an outline of the BNST from a published atlas to confirm the accurate brain region.\n- References for white matter tracts in humans between the BNST and OFC are provided. Is there also evidence for connections between the BNST and the nucleus accumbens and insula?\n- It is unclear how the authors determined that BNST, NAcc, OFC, aInsula were the ROIs of primary interest given that many other regions were significant for the whole-brain analysis. Did they determine with overlap between CUE and TASTE conditions (if so, please show an overlap or conjunction image)? At a threshold? Or were these regions a priori due to prior research?\n- For the self-reported stress levels, please report on the results of the beverage x stress x cue interactions.\n- The authors do a good job modeling the contrasts of interest. Is there enough time between CUE and TASTE to model anticipation of taste for each block?\n- How did the authors identify the threshold for whole brain second levels? Why k = 10?\n- Please provide a rationale for using peaks from the GLM to guide the DCM analysis instead of anatomically defined ROIs, since the introduction is focused on anatomical regions and connections. Please include the details about the size of the spheres used for each brain region in the methods (currently listed in a table).\n- It appears that the authors did not include DCM models for milk vs water in the low stress condition. Please provide a rationale for not investigating connectivity in the low stress condition. It seems that the low stress analysis could provide a baseline measure of effective connectivity and serve as a control for the high stress condition.\n- Were the DCM analyses performed with bilateral regions of interest or separately for each hemisphere? Please specify and include a rationale for the decision.\n- The reviewer understands why the authors excluded participants who rated chocolate milk as less satisfying than water. However, is it possible to perform the DCM model with the participants who rated milk < water, to determine whether their BNST effective connectivity looks different due to their difference in ratings? Second, did any participants have low response rates for ratings and how did the authors handle these low rates?\n- For Figure 4 could the authors please also include the graph for the high stress water cue and high stress water tasted so that the audience can see how those conditions differed from the milk conditions?\n- The authors posit that inhibition of the NAcc via the BNST could promote reward-seeking behaviors to “ameliorate this striatal dampening” that happens during chronic stress. The reviewer thinks this interpretation is over-stated. While the psychosocial stress task used in the study results in dampened NAcc activity, these data are from human neuroimaging studies, where we are unable to confirm the source or underlying biological mechanism of the dampened activity. Can the authors provide a more digestible interpretation of these data that align with the level of measurement we are able to achieve in human neuroimaging studies?\nMinor concerns:\n- Please spell out HS_Cue, etc in figures and text to make it easier for readers to understand the conditions.\n- The reviewer appreciates the authors performing LOO-CV and utilizing the data to also perform predict analyses related to the relation between stress scores and brain effective connectivity. Could the authors reorganize the methods to reflect the organization of the Results section?", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_4": "In this interesting and well-written report, the authors examined effective connectivity between the bed nucleus of the stria terminalis (BNST) and regions involved in reward, goal-directed decision-making, gustatory processing, and interoception during a beverage-tasting task implemented in the context of a brief stressor. Data obtained on a 7T scanner were analyzed using dynamic causal modeling and parametric empirical Bayes, with findings reported as stress-induced inhibitory influence of BNST on cortico-striatal regions (nucleus accumbens, orbitofrontal cortex, anterior insula) during cue exposure and beverage taste. There are several strengths of the work, including a well-reasoned set of hypotheses, a focus on the BNST made possible through inclusion of high-resolution data from 7T, a relatively large sample size of healthy control subjects for this type of task (gustatory processing; n=~48), and rigorous implementation of DCM modeling and PEB approaches. That said, the premise of the study rests on how the study was designed with respect to the stress induction, which unfortunately was methodologically weak. There are also challenges to understanding the specificity of the findings in relation to the stress induction and to rewarding stimuli, due to the analytic approach. Together, these concerns call into question the primary conclusions of the paper. Comments to detail these and related concerns are stated below.\n1) The study design included a standard gustatory (beverage) cue and taste task implemented in the scanner in two runs. Before the first run of the beverage task, a “low-stress” induction was completed (in the scanner), involving 2 minutes of counting backwards from a whole number. Before the second run, a “high-stress” induction was completed (also in the scanner), involving 2 minutes of serial subtractions; errors or responding “too slowly” required starting over. There are several issues with this task and how the task was implemented.\na. The description of the task does not state any of the following critical details for the “high-stress” task: (a) from what number the subtraction started, (b) what numbers were subtracted serially, (c) who gave the feedback, and how this was given (i.e., could the participant see the person giving feedback?), (d) what “too slow” meant in terms of giving feedback on response speed. All of these details are important for interpreting results.\nb. Stressors that are 2 minutes in length are not convincing in terms of actually inducing stress, and the task does not appear to include components of uncontrollability/unpredictability. The negative evaluative threat component is not particularly convincing either. Added to this, objective measures of stress either were not collected (i.e., serum or salivary cortisol) or in the case of electrodermal activity (EDA) did not show any difference between “low-stress” and “high-stress” (p=0.74). A large body of research on robust psychosocial stress tasks, such as the Trier Social Stress Test, the Maastrich Acute Stress Task, and the Montreal Imaging Stress Task are much longer (15-20 minutes), and studies which employ them almost always include an objective measure of HPA-axis reactivity and/or include a well-matched control version which is not subject to order effects (as was the case here, where all subjects first had “low-stress” then “high-stress” stress conditions).\nc. Given the design that did not adequately address the control condition issue, it is nearly impossible to ascertain whether results stemming from the “high-stress” condition were simply a reflection of being in the scanner for as long as they were at that point, having to complete a somewhat challenging task involving cues, water/chocolate milk tastes, swallowing, rinsing, etc.\nd. The references provided as suggesting this task has been validated are insufficiently rigorous as examples for the current study. The first reference (77; Kistenmacher et al) used a very different cross-over design, in which subjects completed 2 visits – one that was control/no stress and one that was stress, and the stress induction was implemented for a full 5 minutes. In addition to these differences, cortisol and ACTH were collected and analyzed to demonstrate a clear stress response in participants during the stress induction visit. The task used in the current study does not match this design.\ne. The other references listed do not provide sufficient detail in their methods to match the current study either, and also did not include objective measures of HPA-axis reactivity, so it is challenging to accept these as “validation” of the task used in the current study.\nf. Although the report does provide analysis of self-reported “stress” ratings, again, given the design, it is not possible to know whether the differences in ratings reported after the “low-stress” induction compared to those reported after the “high-stress” induction simply reflect the influence of demand characteristics, because there is no control arm/visit. This is particularly likely to be the case given the absence of objective, physiological measures showing a stress response (i.e., cortisol/ACTH were not measured, and EDA response showed no difference between conditions).\nUnfortunately, based on the above concerns regarding the design and implementation of the stress induction, the discussion and conclusions are not well supported by the data.\n2) The authors mention that EDA values were “higher after high-stress induction compared to the low-stress induction, though this difference did not reach significance”. This is not only evidence that the stress induction did not actually invoke a physiological response, but is also misleading to state. The mean difference stated was -0.0006 (p=0.74); this is negligible and should not be framed as “higher after high-stress”, even if it didn’t reach significance. Instead, the report of the finding should be revised to make it clear that EDA was not different between conditions, without further elaboration/interpretation.\n3) As it is described, it appears that the DCM model examines beverage cues (chocolate milk + water, or chocolate milk alone) and beverage tastes (again, chocolate milk + water, or chocolate milk alone) specifically in the “high-stress” condition, without comparing to the “low-stress” condition. Therefore, any findings related to the analytic approach prevent inferences related to specificity to stress (especially given concerns described under point 1) and could alternatively have been interpreted as findings specific to completing run 2 of the beverage taste task, not to stress.\n4) Similarly, the model does not appear to contrast chocolate milk to water, in some cases these are combined (models 1 and 3), reducing confidence that the rewarding aspect of the tastants is involved, yet much of the discussion of these findings is framed around reward. This is confusing and perhaps misleading given lack of specificity to the rewarding stimulus (chocolate milk), which was rated as significantly more pleasant than the water stimulus.\n5) A minor point to note is that the authors report excluding 6 subjects whose chocolate milk pleasantness ratings were lower than the water pleasantness ratings, presumably because their ratings would reflect that the chocolate milk was not pleasant, but they do not mention whether data were excluded for any subjects who rated stress as similar across runs. This should be considered because those who rated stress as similar across runs (i.e., those for whom stress ratings did not change) or were lower for the “high-stress” relative to “low-stress” would similarly not fit their assumption that only data for those who found it to be stressful should be included."}, "author_response": "Responses to the Reviewers\nRevised title:\nBed nucleus of the stria terminalis connectivity during high-stress food cue and taste processing\nManuscript: NCOMMS-25-24841\nFor convenience, all comments are provided in boxed format with responses below. Quoted text from the\nmanuscript and supporting information appears in blue. A clickable Table of Contents is provided below to\nfacilitate navigation.\nTable of Contents\nReviewer 1 .................................................................................................................................................2\nReviewer 2 ...............................................................................................................................................12\nReviewer 3 ...............................................................................................................................................25\nReviewer 4 ...............................................................................................................................................26\nThe manuscript by Guerrero-Hreins et al. seeks to address the highly relevant question of whether stress\nmodulates BNST connectivity during taste processing. The study has several strengths, such as using a high\ntemporal resolution 7T sequence that provides a strong resource for dynamic causal modelling, an innovative\napproach based on the experimentally demanding use of an MR-compatible gustometer, and a sufficiently large\nsample. Despite these strengths, there are multiple major issues that require substantial revision before the\nmanuscript can provide the aspired answer. Illustratively, I would expect far more dedication to the right use of\nanatomical regions if this is one of the key contributions of the paper, and I have doubts whether the current\ndefinition is justified.\nResponse:\nWe thank the Reviewer for their thoughtful summary and for acknowledging the strengths of our study. We\nappreciate the reviewer’s emphasis on the importance of anatomical precision, particularly in relation to our\ncentral claims about bed nucleus of the stria terminalis (BNST) connectivity. In our revised manuscript, we have\naddressed these concerns in detail by refining the definition of brain regions and carrying out several new\nanalyses. We trust these revisions have strengthened the manuscript and bring it closer to meeting the\nexpectations outlined.\nMajor issues\n1. Although the introduction is written in an accessible manner, the statement concerning the rationale for\nDCM should be rephrased. I understand that the authors want to argue that this method is better than\n(undirected) FC, but it is only superior to FC if certain assumptions are met (and ideally evaluated). For example,\nit is known that a differential hemodynamic delay can lead to spurious inferences about the directionality of the\neffects. Likewise, differences in tSNR or the seeds can contribute to incorrect inferences. Without additional\ncontrol analyses, it is difficult to take the (strongly worded) inferences for granted.\nResponse:\nWe thank the Reviewer for raising these concerns, which we address—point-by-point—below.\nThe suggestion vis-à-vis DCM’s superiority\nWe do not wish to assert that dynamic causal modelling (DCM) is superior to undirected product-moment\ncorrelation-based approaches in general, but that the former approach is better suited to the task of inferring\ndirected interactions between latent neuronal populations. The rationale has been updated in the manuscript:\nMS, p. 4:\n“However, understanding food cue-related interactions within brain networks has been hampered by the field’s\nreliance on correlational measures (i.e., functional connectivity), rather than methods that can infer effective\nconnectivity—the directed influence that one brain region exerts over another43. Though not superior to simpler\nfunctional connectivity models, dynamic causal modelling (DCM) offers a biophysically grounded framework for\ninferring effective connectivity, and has proven valuable in addressing such questions 44.”\nAdditional control analyses\nRegarding the Reviewer’s other concern, we agree that in silico control analyses are essential for ensuring that\nconnectivity parameters can be recovered in the presence of realistic noise, and we have now included an in\nsilico validation of the main cue- and taste-condition models presented in the manuscript. These analyses\ninvolved (a) estimating subject- and region-specific signal-to-noise ratios (SNRs) for empirical blood-oxygen-\nlevel-dependent (BOLD) time series data, (b) using inferred model parameters (including those that encode\nhemodynamic delays) to simulate BOLD time-series data with region-specific noise (that reflected the estimated\nsimulated data, and evaluating both parameter recovery and the relationship between second-level models\nderived from empirical and synthetic data (see Supporting Information).\nFrom these analyses, we provide a few key takeaways:\n1. Signal quality was consistent and within expected bounds: Estimated SNRs across regions averaged 0.72 ±\n0.01, indicating approximately 39% more background variance than task-related signal.\n2. DCM parameters could be reliably recovered at the subject level: Inversions of noise-matched synthetic data\nproduced maximum a posteriori (MAP) estimates that closely tracked ground-truth values across both cue-\nrelated models (mean r = 0.990 ± 0.012) and taste-related models (mean r = 0.989 ± 0.012).\n3. Group-level inferences were robust to realistic noise and hemodynamic variability: PEB models recovered\nthe pattern of modulatory effects observed in the empirical data, with a root mean squared error (RMSE) of\n0.152 and correlation r = 0.769 for cue-related models, and RMSE of 0.135 and r = 0.994 for taste-related\nmodels.\nBelow, we reproduce Table S7 and Fig. S5. In sum, these results support the validity of our modelling approach\nand the robustness of inferences drawn from empirical data.\nTask modulation Region Mean SNR\nCue BNST 0.72 (0.01)\nCue NAcc 0.73 (0.01)\nCue OFC 0.73 (0.01)\nCue Insula 0.73 (0.01)\nTaste BNST 0.73 (0.01)\nTaste NAcc 0.73 (0.01)\nTaste OFC 0.73 (0.02)\nTaste Insula 0.74 (0.02)\nTable S7. Mean signal-to-noise ratio (SNR) values (± standard\ndeviation) for each region of interest and task modulation. SNR was\ncalculated from the extracted time series used in the DCM analysis and\nindicates consistent signal quality across regions and task conditions.\nRegions include the bed nucleus of the stria terminalis (BNST); nucleus\naccumbens (NAcc); orbitofrontal cortex (OFC); and insula. Cue-related\nDCM models were estimated for 44 participants; taste-related models\nfor 47 participants.\n(PEB-derived) modulatory effects inferred from empirical data (x-axis) and noise-matched synthetic data (y-axis), shown separately\nfor cue-related (a) and taste-related (c) DCMs (cue, n = 44; taste, n = 47). Root mean squared error (RMSE) and Pearson’s correlation\ncoefficient r indicate recovery fidelity. b, d Subject-level parameter recovery assessed by correlating all MAP parameter estimates\n(neuronal and hemodynamic) with ground-truth values across simulated instantiations. Each bar reflects one subject; negative values\nreflect poor identifiability for specific parameter combinations. Recovery is shown for cue (b) and taste (d) conditions, with stronger\nrecovery in the taste dataset overall.\n2. The phrasing of the (clearly) non-significant EDA effects is odd. If the numbers are correct, it is not even\nclose to a directional effect. In my opinion, this should be clearly stated. Since EDA failed to capture\nstress-related physiological effects, the authors may consider using alternatives (e.g., pulse rate as in\n(Kühnel et al., 2022; Erhart et al., 2024)).\nResponse:\nWe thank the Reviewer for the comments, which we address below.\nNonsignificant electrodermal activity result\nWe agree that the characterization of the result is imprecise. We have updated the text to make it clear that this\nresult cannot be considered directionally:\nMS, p. 5:\n“EDA values were not significantly higher after high-stress induction compared to the low-stress induction (t =\n27\n-0.31, p = .76, n = 28, mean difference = -0.006 [-0.05, 0.04]; see Figure S1).”\nSI, Fig. S1:\n“Differences between conditions did not reach statistical significance (t = -0.31, p = .76, n = 28, mean difference\n27\n= -0.006 [-0.05, 0.04]).”\nAlternative physiological measurements\nsuccessive differences (RMSSD) as a measure of heart rate variability (HRV) during the low- and high-stress\nconditions. In accordance with previous research, a window of 20 seconds was used to assess RMSSD\ndifferences between conditions. As expected, mean RMSSD values during the start of the high-stress run were\nlower than during the low stress run t(43) = -2.25, p = 0.029. We provide further details of HRV-related analyses\nin response to a later comment regarding the time-related effects of the stressor.\nSee:\n• Castaldo, R, et al. \"Ultra-short term HRV features as surrogates of short term HRV: A case study on\nmental stress detection in real life.\" BMC medical informatics and decision making 19.1 (2019): 12.\nThis is now detailed in the methods and results section of the manuscript and in the Supporting Information.\nWe thank the reviewer again for their helpful suggestion, which has contributed valuable physiological evidence to\nsupport our stress induction procedure.\nMS, p. 5:\n“Mean RMSSD HRV values during the beginning of the high-stress run were lower than during the beginning of\nthe low stress run (t43 = -2.25, p = 0.029), indicating elevated physiological stress response.”\nMS, p. 16:\n“The root mean square of the successive differences (RMSSD) in pulse signals extracted using the PhysIO\ntoolbox were used as a measure of HRV during the low- and high-stress conditions. In accordance with previous\nresearch102,103, a window of 20 seconds was used to assess RMSSD differences at the beginning of each run.\nTo evaluate the sustained effects of stress induction, we analysed the complete time course of RMSSD during\nthe beverage task using a sliding-window approach (10-second windows with 2-second steps). RMSSD values\nwere extracted continuously across the task period and subsequently segmented into early, middle, and late\nphases, defined as the first, middle, and final 20 seconds of the task, respectively. These values were entered\ninto a within-participants repeated-measures ANOVA with Condition (low-stress, high-stress) and Time (early,\nmiddle, late) as factors. Post hoc comparisons between time points within each condition were conducted using\npaired-sample t-tests to assess directional trends over the course of the task.”\n3. In general, the use of causal language for functions of brain regions (“OFC is responsible for”) is not\njustified unless there is very clear evidence from lesion studies. I strongly suggest rephrasing these\nsections to tone down these statements.\nResponse:\nWe fully agree with the reviewer’s comment and have considerably changed the language used in the entire\nmanuscript to avoid unwarranted causal language. For example, the sentence regarding OFC function in the\ndiscussion has been changed and both the introduction and discussion have undergone substantial changes in\ntheir tone and content. All changes to the manuscript text have been highlighted.\nMS, p. 7:\n“OFC activity has been associated with the encoding of reward predictions necessary for value-based goal-\ndirected action.”\n4. The labels for the seed regions are highly misleading in two cases. First, the selected seed is not in the\nanterior insula, but close to the dorsal mid insula, likely including the primary ingestive cortex. This is\nplausible given the design of the task. Still, it weakens the argument that direct structural connections to\nthe BNST exist (they are mostly to the anterior part of the insula, which is also involved in taste\nprocessing, but likely more concerning expectations). Second, the BNST highlighted in the figures does\nnot correspond well to the BNST. The reference to the Pauli atlas is also misleading in that sense because\nthey only provide a mask for the extended amygdala. This includes the BNST, but also other parts, and\nit seems as if the authors have used the entire mask as search volume, leading to a situation where the\ngroup-level peaks are located well outside the BNST. Also, the provided coordinate numbers do not\ncorrespond to the figure, so I guess a minus sign is missing for the z coordinate (otherwise, they would\nbe placed even further away from the targets). In my opinion, this is a major problem since the paper is\nfocused on the function of the BNST while the current analysis applies more to the extended amygdala,\nwhich is much more difficult to connect to the preclinical work.\nResponse:\nWe address each point raised by the Reviewer below:\nInsula seed label clarification\nWe thank the Reviewer for this observation. We agree that the label anterior insula is not precise, and we have\nrevised all labels in the manuscript and figures to refer to this region as the insula rather than anterior insula.\nFurther perusal of the Flook et al.’s ‘BNST-insula structural connectivity in humans’ publication shows that our\ninsula VOI falls within the subregion which the authors labelled as the middle short gyrus. This subdivision,\ncompared to posterior subdivisions, was found to have significant connectivity with the BNST. We have changed\nthe text in the manuscript to properly describe this:\nMS, p. 8:\n“Structural connectivity between the BNST and insula has been recently confirmed in humans using probabilistic\ntractography31, with the mid-insula subdivision used as our DCM VOI found to have significant connectivity with\nthe BNST, albeit to a lesser extent than anterior insula subdivisions.”\nBNST mask and peak location\nWe thank the Reviewer for their observations. We agree that the Pauli extended amygdala atlas encompasses\na larger area than the BNST alone. To address this issue, we have rerun our DCM analyses using a new BNST\nmask, created by overlaying a manually segmented BNST anatomical mask from Theiss et al., (2017) onto\nPauli’s extended amygdala mask (Pauli et al., 2018). Both masks were binarised and combined to create a new\nconsensus mask containing only voxels common to the original masks. We have also updated the BNST z-\ncoordinate number in Table 1, which was incorrectly missing a minus sign as Reviewer 2 also pointed out.\nSee:\n• Theiss, J. D., Ridgewell, C., McHugo, M., Heckers, S., & Blackford, J. U. (2017). Manual segmentation\nof the human bed nucleus of the stria terminalis using 3T MRI. NeuroImage, 146, 288–292.\nOur obtained results using this new BNST mask remain largely consistent. For example, evidence of BNST\ninhibition on the other regions in our model space during high-stress cue presentation was found using this\nconsensus BNST mask. However, not all the originally observed high-stress taste effects remained significant.\nThese changes are now reflected throughout the manuscript.\nThe updated Table 1 and Figure 2 have been reproduced below:\nTable 1\nRegion MNI coordinates\nx y z\nBNST 22 -8 -11\nNAcc 10 19 -8\nOFC 43 42 -19\nINS 37 5 10\nBNST, bed nucleus of the stria terminalis; INS, insula; MNI, Montreal Neurological Institute;\nNAcc, nucleus accumbens; OFC, orbitofrontal cortex. Time-series extraction for the BNST\nwas run using a consensus mask, created by overlaying the manually segmented BNST\nmask from Theiss and colleagues95, onto Pauli and colleagues’ extended amygdala mask96.\nBoth masks were binarized and combined to create a new mask containing only voxels\ncommon to both original masks. An inclusive mask was created for the NAcc using the 1-\nmm3 resolution atlas from the MNI standardized automated anatomical labelling 3\ntoolbox (AAL3,98). The first eigenvariate across all activated voxels within a 4 mm radial\nsphere around the subject-specific maxima (p < .05, uncorrected), but no more than 8 mm from the\ngroup maxima, were extracted.\nFigure 2. Behavioural ratings and group-level whole-brain general linear model activation\nduring food cue and taste processing. (a) Mean self-reported hunger and stress ratings during\nlow-stress and high-stress inductions, paired t-tests, *p < .05, ****p < .0001 (low stress vs high\nstress). Error bars show standard error of the mean. (b) Mean self-reported chocolate milk and\n7\nwater pleasantness ratings after low- and high-stress inductions, 2-way ANOVA (stress x\nbeverage), main effect of beverage, ****p < .0001, post hoc multiple comparisons, paired t-tests,\n*p < .05, **p < .01 (low-stress vs high-stress). (c) Time course showing mean self-reported\nchocolate milk and water pleasantness ratings during each beverage cue exposure, 3-way\nANOVA (stress x beverage x cue), post hoc multiple comparisons unpaired t-tests at each cue\nexposure, ***p < .001, ****p < .0001 (chocolate milk vs water during high stress); ^^^p < .001;\n^^^^p < .0001 (chocolate milk vs water during low stress); n = 48 for all. Error bars show standard\nerror of the mean. (d) Group-level GLM results for cue presentation (CUE), n = 44. (e) Group-\nlevel GLM results for taste receipt (TASTE), n = 47. Group-level activation t-maps were overlaid\non a standard MNI space template. The t-values are represented with the colour bar. Significant\nactivation was thresholded at a voxel-wise pFDR < .05, KE ≥ 10. Voxels of interest are highlighted\nwith yellow-dotted circles, including the OFC, INS, BNST and NAcc. An inclusive binary mask\nwas used during the NAcc time series extraction from the MNI-standardized automated\nanatomical labelling 3 (AAL3, 98) toolbox, shown in red, and BNST time series extraction (see\nMethods) shown in yellow. BNST, bed nucleus of the stria terminalis; INS, insula; NAcc, nucleus\naccumbens; OFC, orbitofrontal cortex; HS, high stress; LS, low stress.\n5. It is not clear to me what the evidence is that BNST input would inhibit goal-directed behavior as there\nwas no association with behavior that would support this (should be rephrased). In general, the authors\ndo not make much use of the behavioral ratings to strengthen their behavioral claims in the discussion. I\nthink they should provide additional (exploratory) analyses with changes in stress and liking or edit the\n(far too strongly worded) language in the discussion.\nResponse:\nWe agree with the reviewer’s comment and have removed text suggesting that the BNST inhibits goal-directed\nbehaviour within the context of this task. As recommended by the reviewer, we have rewritten this section to\nunderscore the behavioural findings that support that this pathway is linked to stress sensitivity and cue\nprocessing – factors directly related to the task itself. This section in the discussion now reads as follows:\nMS, p. 7:\n“OFC activity has previously been associated with the encoding of reward predictions necessary for valence-\nbased action 48–51 and OFC-mediated shifts in behaviour under stressful conditions33 likely involves homeostatic\ninput from sensory regions, such as the BNST14, as extended amygdala lesions have been found to result in\ninaccurate cue-outcome predictions 52. We also found that the strength of BNST downregulation of the OFC\npredicted increased self-reported stress sensitivity during the task. This indicates that BNST-mediated\nmodulation of cortical regions may be sensitive to experiential stress levels, which could, in turn, viably promote\nreinforced behaviours, such as food approach53. Further investigation on the functional role of this pathway using\nparadigms that compare goal-directed versus habitual food approach under stress is warranted.”\nAs suggested, we have also significantly overhauled the tone of the discussion to tone down the inferences that\nare made from our obtained results. These changes are highlighted in yellow in the revised manuscript.\n6. The discussion should be rebalanced to better reflect which findings speak against the (posited)\nspecificity. For example, isn’t it surprising in light of the hypothesis that beverage taste receipt under\nstress elicited similar inhibitory effects and was not specific to the milkshake? How would the authors\nrule out that they are reporting a less functionally distinct process here? Water is essentially also a control\ncondition for everything that is simply driven by the timing of the task elements, and I do not think this\ninterpretation was given enough consideration in the report. This concern is even more important given\nthe lack of randomization of the block order.\nResponse:\nWe agree that the observation of similar inhibitory effects during taste receipt for both chocolate milk and water\nunder high stress raises important questions about specificity. Although our original hypothesis focused on\nstress-mediated modulation of palatable food processing, our findings during taste receipt suggest a broader\nstress-induced effect on BNST-insula connectivity that is not exclusive to reward value. We now emphasize that\nwater stimuli under high stress is not solely a neutral control in this context as it is contrasted to the low-stress\nwater stimuli in our model. Importantly, the cue and taste conditions involve sensory stimulation, oral motor\nresponse, and task expectancy, all of which may engage insula-BNST circuitry. The insula's role in integrating\ninteroceptive and sensory signals could explain its involvement irrespective of reward content, particularly under\nstress when interoceptive sensitivity is altered. As suggested, we have changed the discussion to highlight these\npoints.\nMS, p. 9:\n\"Both the cue and taste phases of the task involved sensory input, oral motor engagement, and expectancy—all\nfeatures known to recruit insula-BNST pathways66,67. Considering the insula’s role in integrating interoceptive\nand sensory signals, its engagement could plausibly reflect general stress-related modulation rather than reward-\nspecific processing, particularly under conditions of altered interoceptive sensitivity induced by acute stress.”\nRegarding block order, we used a fixed low-stress and high-stress sequence to minimize carryover effects,\nconsistent with prior work. While this limits counterbalancing, it avoids contamination of the low-stress baseline\nwith residual stress effects. We have added text to the limitations section stating that future studies should\nimplement washout periods or counterbalanced designs to further isolate stress-specific contributions.\nMS, p. 9:\n“Future studies should aim to implement washout periods or counterbalanced designs to better isolate stress-\nspecific contributions.”\n7. Please restructure the information so that the paper can be read linearly. For example, it is only later\nrevealed why the authors excluded 3 participants for the analysis of the cue responses? Were the\nexclusion criteria prespecified, and would the TASTE DCM with the same 45 participants give the same\nresults? Wouldn’t this approach inflate the model estimates if people without a clear activation in VOIs\nare selectively excluded (was it due to the BNST or various regions)? Couldn’t the subject-specific peak\nfrom the TASTE analysis be used instead?\nResponse:\nWe thank the Reviewer for this helpful observation. We have revised this information to be included earlier in the\nmanuscript so that it can be read linearly. Please note that the number of excluded participants has slightly\nincreased from 3 to 4 participants as a result of the new consensus BNST mask.\nMS, p. 6:\n“Four participants were excluded from the ‘CUE’ DCM analyses (n = 44) and one from the ‘TASTE’ DCM\nanalyses (n = 47) as they did not survive the BNST VOI time-series extraction.”\nunderstandable. The exclusion criteria for VOIs adhered to the guidelines suggested by the developers of DCM.\nSee:\n• Zeidman, Peter, et al. \"A guide to group effective connectivity analysis, part 1: First level analysis with\nDCM for fMRI.\" Neuroimage 200 (2019): 174-190.\nAs outlined in Appendix 1 of the publication, participants were not originally excluded for lack of suprathreshold\nresponses. Instead, when no voxel met the initial threshold, it was progressively lowered until a peak could be\nidentified up to a threshold of 0.15. This approach avoids selection bias while maintaining anatomical specificity.\nMoreover, DCM’s Bayesian framework penalizes model complexity and incorporates shrinkage priors, so\nparameter estimates remain conservative in the absence of strong signal. However, these excluded participants\ndid not produce any significant voxels in the BNST at a very liberal p value threshold of 0.15 during CUE,\nindicating that reliable signal from the BNST cannot be acquired for these individuals during these specific\nconditions. As such, these participants cannot be included in the analyses as they do not produce a viable time\nseries. Likewise, using peaks from another condition would be unprincipled within the DCM framework, as it\nseparates the link between DCM model specification and the experimental effect under study.\n8. Perhaps I missed it, but I do not see any attempt to account for the potential time-related effect of the\npreceding stressor. I understand that it is difficult to dissociate linear trends from signal drift in the MR\nsignal, but the beverage runs are quite long compared to the preceding stress exposure. Was the\navailable physiological data (e.g., changes in pulse rate) explored to evaluate how long the after-effects\nof stress last?\nResponse:\nThe Reviewer raises an important point about how pulse rate trace can be used to evaluate the effectiveness of\nour stress induction task. As mentioned above—in response to the second comment—we have evaluated\nRMSSD pulse signals as a measure of HRV during the low- and high-stress conditions.\nMean RMSSD values during the start of the high-stress run were lower than during the low-stress run, t = -\n(43)\n2.25, p = 0.029. To directly address the Reviewer’s suggestion to evaluate the sustained after-effects of high\nstress induction, we analyzed the complete time course of RMSSD using a sliding-window approach (10-second\nwindows stepped every 2 seconds) across the beverage task. We then segmented the data into early, middle,\nand late task periods (first, middle, and final 20 seconds) to evaluate trends across time. The within-participants\nrepeated-measures ANOVA showed a significant Condition × Time interaction F = 3.20, p = 0.045, indicating\n(2, 86)\nthat RMSSD evolved differently over time between the low- and high-stress conditions. Additional paired t-test\ncomparisons showed that RMSSD increased significantly more from early to middle in the LS condition\ncompared to the high-stress condition, t = 2.25, p = 0.029. Similarly, from early to late, RMSSD increased in\n(43)\nthe low-stress condition but decreased in the high-stress condition, t = 2.03, p = 0.048.\n(43)\nThese findings suggest that the physiological impact of the stressor persisted during beverage task, with limited\nevidence of recovery. This pattern supports that the after-effects of the stress manipulation were not transient\nand may have influenced participants’ neural response throughout the task. These analyses are now reported\nin the manuscript:\nMS, p. 5:\n“Relatedly, the complete time course of RMSSD across both task runs was analysed using a sliding-window\napproach (10-second windows stepped every 2 seconds). Data were segmented into early, middle, and late task\nperiods (first, middle, and final 20 seconds) to evaluate trends across time. A within-participants repeated-\nmeasures ANOVA found a significant Condition × Time interaction F = 3.20, p = 0.045, indicating that\n(2, 86)\nRMSSD evolved differently over time between the low- and high-stress conditions. Additional paired t-test\ncomparisons showed that RMSSD increased significantly more from early to middle in the low-stress condition\ncompared to the high-stress condition, t = 2.25, p = 0.029. Likewise, RMSSD increased in the low-stress\n(43)\ncondition but decreased in the high-stress condition when comparing early to late windows, t = 2.03, p =\n(43)\n0.048.”\n9. The sequence is optimized for acquisition speed, and this may improve the effective connectivity\nestimates. However, in my experience, the use of GRAPPA and a multiband factor >4 can substantially\nreduce the SNR in deeper brain regions. Did the authors calculate and evaluate tSNR in their target\nregions, since differences in tSNR (e.g., between the insula and the BNST) could also affect the model\nestimates of effective connectivity (since BNST signal will likely fluctuate much more).\nResponse:\nWe thank the Reviewer for this important point. We agree that high multiband acceleration factors reduce\ntemporal SNR. As outlined above (in response to the first issue raised), we estimated region-specific SNR directly\nfrom empirical BOLD data (including both cortical and subcortical regions of interest). These empirical estimates\nwere then used to inform our in silico validation analyses. As reported in the Supporting Information (Figures S4\n& S5), our parameter recovery results demonstrate that effective connectivity estimates remained robust.\nRemarks on code availability:\nI have not reviewed the code yet because I think the definition of the VOIs is a crucial aspect that needs to be\naddressed first (and will likely lead to different results).\nResponse:\nAs outlined in our response above, we have now clarified the VOI definitions, updated region labels where\nappropriate, and revised the interpretation of affected results. The updated results are now available on the\nstudy’s github repository: https://github.com/evaghreins/Gustometer_Stress_BNST_DCM\nWe welcome any further feedback on the code now that these clarifications have been made.\nThe authors present a fascinating paper related to the role of the bed nucleus of the stria terminalis (BNST) in\nstress-related feeding. The authors used a form of psychosocial stress to induce stress in participants, along\nwith CUE and TASTE conditions related to the receipt of water or chocolate milk, respectively. The authors also\nused dynamic causal modeling to better determine direction of inhibition and excitation within a network\ncomprised of the BNST, nucleus accumbens (NAcc), orbitofrontal cortex (OFC), and anterior insula (aInsula).\nThe strengths of the paper are rationale based on good evidence from rodent models, strong experimental task,\nand generally excellent methods including covarying for physiological signal and leave one out analyses. While\nthis is a great paper, there are several major and minor concerns related to the methodology and discussion of\nthe data that reduce enthusiasm.\nResponse:\nWe thank the Reviewer for their thoughtful and encouraging summary of our work, and for recognizing the\nstrengths of our experimental design and modeling approach. We appreciate the constructive concerns raised\nand have revised the manuscript to address each point in detail—clarifying our methods and expanding on key\naspects of the discussion.\nMajor concerns:\n- The major concern is the definition of the brain regions used in the analyses. The area of activation listed\nas BNST in Figure 1 appears to be dorsal amygdala and not BNST. For the coordinates shown in Table 1, maybe\nthere is an error for the BNST peak (should be z = -11?); the peak coordinates for the NAcc appears to be in the\nvery most anterior portion, very close to white matter; the peak coordinates for the anterior insula appear to be\nin the dorsal middle or posterior insula (see 10.1093/scan/nss066). It would be helpful to overlay the areas used\nin the DCM analysis on Figure 1. For the enlarged parts of Figure 1, please show an outline of the BNST from a\npublished atlas to confirm the accurate brain region.\nResponse:\nWe address each point raised by the Reviewer below:\nBNST mask and peak location\nWe thank the Reviewer for their observation. We agree that the mask used, taken from the Pauli extended\namygdala atlas, encompasses a larger area than the BNST alone. To address this issue, we have rerun our\nDCM analyses using a new BNST mask, created by overlaying a manually segmented BNST anatomical mask\nfrom Theiss et al., (2017) onto Pauli’s extended amygdala mask (Pauli et al., 2018). The masks were binarised\nand combined, resulting in a new limited mask containing only voxels common to both original masks. We have\nalso updated the BNST z-coordinate number in Table 1 to -11, which was missing a minus sign as the Reviewer\ncorrectly pointed out.\nAs suggested by the Reviewer, we have also included a figure showing our new BNST mask overlaid with the\nTheiss BNST mask (2017) on an MNI template, provided below.\nThe updated Table 1 and Figure 2 have been reproduced below:\nGroup peak coordinates for DCM volumes-of-interest\nRegion MNI coordinates\nx y z\nBNST 22 -8 -11\nNAcc 10 19 -8\nOFC 43 42 -19\nINS 37 5 10\nBNST, bed nucleus of the stria terminalis; INS, insula; MNI, Montreal\nNeurological Institute; NAcc, nucleus accumbens; OFC, orbitofrontal cortex.\nTime-series extraction for the BNST was run using a consensus mask, created\nby overlaying the manually segmented BNST mask from Theiss and\ncolleagues95, onto Pauli and colleagues’ extended amygdala mask96. Both\nmasks were binarized and combined to create a new mask containing only\nvoxels common to both original masks. An inclusive mask was created for the\nNAcc using the 1-mm3 resolution atlas from the MNI standardized automated\nanatomical labelling 3 toolbox (AAL3,98). The first eigenvariate across all\nactivated voxels within a 4 mm radial sphere around the subject-specific\nmaxima (p < .05, uncorrected), but no more than 8 mm from the group maxima,\nwere extracted.\nfood cue and taste processing. (a) Mean self-reported hunger and stress ratings during low-stress\nand high-stress inductions, paired t-tests, *p < .05, ****p < .0001 (low stress vs high stress). Error bars\nshow standard error of the mean. (b) Mean self-reported chocolate milk and water pleasantness ratings\nafter low- and high-stress inductions, 2-way ANOVA (stress x beverage), main effect of beverage, ****p\n< .0001, post hoc multiple comparisons, paired t-tests, *p < .05, **p < .01 (low-stress vs high-stress).\n(c) Time course showing mean self-reported chocolate milk and water pleasantness ratings during each\nbeverage cue exposure, 3-way ANOVA (stress x beverage x cue), post hoc multiple comparisons\nunpaired t-tests at each cue exposure, ***p < .001, ****p < .0001 (chocolate milk vs water during high\nstress); ^^^p < .001; ^^^^p < .0001 (chocolate milk vs water during low stress); n = 48 for all. Error bars\nshow standard error of the mean. (d) Group-level GLM results for cue presentation (CUE), n = 44. (e)\nGroup-level GLM results for taste receipt (TASTE), n = 47. Group-level activation t-maps were overlaid\non a standard MNI space template. T-values are represented with the colour bar. Significant\nactivation was thresholded at a voxel-wise pFDR < .05, KE ≥ 10. Voxels of interest are highlighted with yellow-\ndotted circles, including the OFC, INS, BNST and NAcc. An inclusive binary mask was used during the\nNAcc time series extraction from the MNI-standardized automated anatomical labelling 3 (AAL3, 98)\ntoolbox, shown in red, and BNST time series extraction (see methods) shown in yellow. BNST, bed\nnucleus of the stria terminalis; INS, insula; NAcc, nucleus accumbens; OFC, orbitofrontal cortex; HS,\nhigh stress; LS, low stress.\nConsensus BNST mask (Green) overlaid with the Theiss BNST (2017) mask (Red) on an MNI template\n14\nWe agree that the seed used encompasses the dorsal mid-insula as opposed to the anterior insula. We have\nrevised all labels in the manuscript and figures to refer to this region as the insula rather than anterior insula. We\nhave included a figure with the high-stress CUE and TASTE peak coordinates of all participants for all four\nregions.\nFigure S3. Individual participant peak coordinates during (a) high stress cue and (b) high stress taste across the\nBNST, NAcc, INS, and OFC. Each sphere represents the peak voxel for one participant and is coloured\naccording to the region of interest. Peaks are displayed in MNI space on a standard brain surface using BrainNet\nViewer. Overlapping points reflect spatial consistency across participants. BNST = bed nucleus of stria\nterminalis; INS = insula; NAcc = nucleus accumbens; OFC = orbitofrontal cortex.\nOur preprocessing pipeline included aCompCor to model and remove noise from white matter and CSF, reducing\nthe risk of signal contamination in nearby regions. This, combined with DARTEL-based segmentation, ensured\nthat the NAcc mask is restricted to grey matter and excludes adjacent white matter or CSF.\n- References for white matter tracts in humans between the BNST and OFC are provided. Is there also\nevidence for connections between the BNST and the nucleus accumbens and insula?\nWe thank the reviewer for their comment. There is published evidence for human white matter track connections\nbetween the BNST and insula (Flook et al., 2020), and between the BNST and the nucleus accumbens (Kruger\net al., 2015) using diffusion-weighted MRI and probabilistic tractography.\nSee:\n• Flook, E. A. et al. BNST-insula structural connectivity in humans. Neuroimage 210, 116555 (2020).\n• Krüger, O., Shiozawa, T., Kreifelts, B., Scheffler, K. & Ethofer, T. Three distinct fiber pathways of the bed\nnucleus of the stria terminalis to the amygdala and prefrontal cortex. Cortex 66, 60–68 (2015).\nWe have made this clearer in the manuscript:\nMS, p. 3:\n“Evidence from diffusion-weighted MRI and post-mortem studies have confirmed bidirectional white matter tracts\nbetween the BNST, NAcc, insula and OFC in humans28–31, such that these regions form an anatomical network.”\nMS, p. 9:\n“Structural connectivity between the BNST and insula has been recently confirmed in humans using probabilistic\ntractography31...”\n- It is unclear how the authors determined that BNST, NAcc, OFC, aInsula were the ROIs of primary\ninterest given that many other regions were significant for the whole-brain analysis. Did they determine\nwith overlap between CUE and TASTE conditions (if so, please show an overlap or conjunction image)?\nAt a threshold? Or were these regions a priori due to prior research?\nResponse:\nWe thank the Reviewer for their observation. The BNST was an a priori region of interest based on the emerging\nanimal literature delineating a functional role of the BNST in stress, reward and feeding (e.g., Daneil, et al., 2015;\nde Araujo Salgado, 2023). The NAcc, OFC and insula ROIs were chosen based off both the CUE and TASTE\nGLM results, which demonstrated that activity in these regions was modulated by our task conditions. Such an\napproach based upon including regions from obtained GLM results is suggested in DCM guidelines (Zeidman et\nal., 2019). These regions were also included in our DCM model space based on their known functional relevance\nwith regards to stress, reward and feeding.\nAs suggested, we have included a labelled overlapped image of the CUE and TASTE group-level activation\nshowing all four ROIs at a threshold of p < 0.05 and included it in the Supporting Information:\nFDR\nSee:\n• Daniel, S. E. & Rainnie, D. G. Stress Modulation of Opposing Circuits in the Bed Nucleus of the Stria\nTerminalis. Neuropsychopharmacol 41, 103–125 (2016).\n• De Araujo Salgado, I. et al. Toggling between food-seeking and self-preservation behaviors via\nhypothalamic response networks. Neuron 111, 2899-2917.e6 (2023).\n• Zeidman, P. et al. A guide to group effective connectivity analysis, part 1: First level analysis with DCM\nfor fMRI. NeuroImage 200, 174–190 (2019).\nregions of interest (ROIs). Overlaid statistical maps showing thresholded group-level activation (FDR-\ncorrected, p < 0.05) for high-stress cue (shown in orange) and high-stress taste (shown in blue)\nconditions across the bed nucleus of stria terminalis (BNST), insula (INS), nucleus accumbens (NAcc)\nand orbitofrontal cortex (OFC). ROIs are shown on a standard brain surface in MNI space. Colour\noverlays indicate spatial extent and overlap of activation between tasks within each ROI. BNST Taste\nactivation outlined in white, NAcc Taste activation outlined in red. BNST = bed nucleus of stria\nterminalis; INS = insula; NAcc = nucleus accumbens; OFC = orbitofrontal cortex.\n- For the self-reported stress levels, please report on the results of the beverage x stress x cue interactions.\nResponse:\nWe thank the Reviewer for your suggestion. The results of the beverage x stress x cue interactions have been\nadded to the main body of the results section, in addition to the legend of Figure 2. This has been reproduced\nbelow.\nMS, p. 5:\nSelf-reported stress levels were higher after the high-stress induction (mean 5.04, standard deviation [SD] 1.58)\ncompared to the low-stress induction (mean 2.27 [SD 1.48],t = 12.41, p < .0001, n = 48; Figure 2a). There\n(47)\nwas a main effect of stress (F = 15.01, p = .0002, n = 48) and beverage (F = 38.91, p < .0001, n = 48)\n(1,94) (1,94)\non pleasantness ratings: chocolate milk was rated more pleasant than water after both high-stress (mean 5.26\n18\n[SD 1.23] vs mean 3.88 [SD 1.26]) and low-stress induction (mean 5.61 [SD 1.12] vs mean 4.29 [SD 1.04], p <\n.0001 for both; Figure 2b, Supporting Information). There was a main effect of cue exposure (F = 10.17, p\n(5,940)\n< .0001), an interaction between cue exposure and stress (F = 2.64, p = .02), and no significant interaction\n(5,940)\nbetween stress, beverage type and cue exposure (F = 1.007, p = 0.41, Figure 2c). Post hoc multiple\n(5,940)\ncomparisons unpaired t-tests revealed chocolate milk was rated more pleasant than water at each cue exposure\nduring high stress (p < .001 at cue exposure 1, 2, 3 and 5; p < .0001 at cue exposure 4 and 6, n = 48 for all) and\nlow stress (p < .001 at cue exposure 4, p < .0001 at cue exposure 2, 3, 5, and 6, n = 48 for all). These results\nconfirm that (1) the high-stress manipulation was successful in increasing participants’ self-reported stress levels,\nand (2) participants rated the rewarding beverage (chocolate milk) as more pleasant than the neutral beverage\n(water) throughout the beverage task.\n- The authors do a good job modeling the contrasts of interest. Is there enough time between CUE and\nTASTE to model anticipation of taste for each block?\nResponse:\nWe thank the Reviewer for this thoughtful suggestion. Modelling anticipation separately from cue presentation\nand taste delivery is of interest; however, in our current design, the 5-second window for the CUE and TASTE\nevents is insufficient to reliably isolate the anticipatory neural responses. Our task design and the nature of the\nhemodynamic response prevent us from adequately disentangling closely spaced events without introducing the\nrisk of collinearity or model instability. We agree that future studies could extend the anticipation window to model\nthis component.\n- How did the authors identify the threshold for whole brain second levels? Why k = 10?\nResponse:\nWe used voxel-wise FDR p < 0.05 correction to determine statistical significance. This approach is more\nappropriate than FWE for 7T fMRI due to the high number of voxels (resels) that reduced FWHM smoothing and\nhigher resolution deliver. FWE, for example, has a very high risk of type II error when reduced smoothing is\napplied, whereas FDR adequately controls for type I error without sacrificing sensitivity.\nSee:\n• Viessmann, O., & Polimeni, J. R. (2021). High-resolution fMRI at 7 Tesla: challenges, promises and\nrecent developments for individual-focused fMRI studies. Current opinion in behavioral sciences, 40, 96–\n104.\n• Tierney, T.M., Clark, C.A. and Carmichael, D.W., 2016. Is Bonferroni correction more sensitive than\nrandom field theory for most fMRI studies?. arXiv:1607.08205\nThe k = 10 cluster extent threshold was used for visualization and table reporting purposes only, and this is now\nspecified in the Methods section:\nMS, p. 13:\nA false discovery rate (FDR) corrected statistical threshold was applied (PFDR<.05) to all whole-brain GLM\nanalyses with a 10-voxel cluster-extent threshold (KE ≥10 voxels). A KE ≥10 extent threshold was used for\nvisualization and reporting purposes.\n- Please provide a rationale for using peaks from the GLM to guide the DCM analysis instead of\nanatomically defined ROIs, since the introduction is focused on anatomical regions and connections.\n(currently listed in a table).\nResponse:\nWe thank the Reviewer for this comment as we agree that it would be useful for readers to understand why these\npeaks were selected and how the ROI masks were created.\nAs outlined in DCM analysis guidelines from Zeidman and colleagues, GLM-derived peak activations should be\nused to define regions of interest (ROIs) for time series extraction. This approach aims to ensure that the\nextracted signal reflects task-relevant neural activity. All the regions included in the DCM network model space\nwere significant in the GLM contrasts and our hypotheses were also anatomically motivated in that we only\nincluded regions that have been demonstrated to show structural connectivity to the BNST. We have now\nclarified this rationale in the manuscript:\nMS, p. 6:\n“As suggested by published guidelines for DCM45,46, the selection of VOIs for analysis was informed by GLM\nactivation results to ensure that VOIs contained task-relevant activity.”\nThe following information has been added to the methods section regarding how masks were created:\nMS, p. 14:\n“For the insula and OFC, the first eigenvariate across all activated voxels within a 4 mm radial sphere around\nthe subject-specific maxima (p < .05, uncorrected), but no more than 8 mm from the group maxima, were\nextracted. Time-series extraction for the BNST was run using a consensus mask, created by overlaying the\nmanually segmented BNST mask from Theiss and colleagues95, onto Pauli and colleagues’ extended amygdala\nmask96. Both masks were binarized and combined to create a new mask containing only voxels common to the\noriginal masks. The NAcc mask was generated using the 1-mm3 resolution atlas from automated anatomical\nlabelling 3 (AAL397) toolbox. The OFC and insula were modelled unilaterally, using the strongest task-evoked\npeaks from the GLM, in line with DCM guidelines prioritizing dominant hemispheric activation. The BNST and\nNAcc were modelled bilaterally given their small size, proximity to the midline, and known bilateral function.”\n- It appears that the authors did not include DCM models for milk vs water in the low stress condition.\nPlease provide a rationale for not investigating connectivity in the low stress condition. It seems that the\nlow stress analysis could provide a baseline measure of effective connectivity and serve as a control for\nthe high stress condition.\nResponse:\nWe thank the Reviewer for your comment.\nAlthough the primary focus of our analysis is on effective connectivity during the high-stress condition in\ncomparison to the low-stress condition, we acknowledge the value to readers of having connectivity values from\nthe low-stress run. We have run separate DCM models for Cue and Taste effective connectivity during the low-\nstress condition and included the results in the Supporting Information. Although this analysis falls outside the\nmain scope of the paper as our choice to focus on the BNST was based upon its role in stress reactivity, we\nagree it provides useful context and helps to motivate future work investigating BNST network connectivity.\nPlease find a summary table of the results reproduced below:\nTable S3. Modulatory influence of low-stress cues and taste\nLow Stress Cue INS --> BNST 0.79 1.00\nBNST --> NAcc 0.40 1.00\nBNST --> OFC 0.64 1.00\nBNST --> INS 0.34 0.98\nLow Stress Chocolate Milk Cue INS --> BNST 0.51 0.82\nOFC --> BNST 0.45 0.75\nBNST --> NAcc 0.30 0.97\nBNST --> OFC 0.41 0.99\nLow Stress Taste INS --> BNST 0.63 1.00\nBNST --> NAcc 0.32 1.00\nBNST --> OFC 0.17 0.74\nBNST --> INS 0.23 0.91\nLow Stress Chocolate Milk Taste OFC --> BNST 0.53 0.92\nBNST --> NAcc 0.30 0.99\nBNST --> INS 0.22 0.77\nBNST, bed nucleus of the stria terminalis; Hz, hertz; INS, insula; MAP, maximum a posteriori estimate; NAcc,\nnucleus accumbens; OFC, orbitofrontal cortex. Pp, posterior probability. The Pp threshold cutoff was set at 0.75\n(moderately strong evidence). Continuous-time modulation is framed in terms of rates (measured in Hz),\nwhereby a task (e.g. Low Stress Cues) modulates the degree to which activity of an efferent node affects the\ncontinuous rate at which the state of afferent node changes.\n- Were the DCM analyses performed with bilateral regions of interest or separately for each hemisphere?\nPlease specify and include a rationale for the decision.\nResponse:\nThe OFC and insula were modeled unilaterally, based on strongest task-evoked peaks from the GLM, to focus\non dominant hemispheric activation as suggested in guidelines by the developers of DCM (Zeidman et al, 2018).\nIn contrast, the BNST and NAcc were modeled as a single bilateral ROI due to their small size, midline-adjacent\nlocation, and well-established bilateral functional organization. As suggested by the reviewer, we have provided\nour rationale for this approach in the manuscript:\nMS, p. 14\n“The OFC and insula were modeled unilaterally, using the strongest task-evoked peaks from the GLM, in line\nwith DCM guidelines prioritizing dominant hemispheric activation. The BNST and NAcc were modelled bilaterally\ngiven their small size, proximity to the midline, and known bilateral function.”\n- The reviewer understands why the authors excluded participants who rated chocolate milk as less\nsatisfying than water. However, is it possible to perform the DCM model with the participants who rated\nmilk < water, to determine whether their BNST effective connectivity looks different due to their difference\nin ratings? Second, did any participants have low response rates for ratings and how did the authors\nhandle these low rates?\nResponse:\nWe thank the Reviewer for your comment.\nAs suggested, the DCM models for CUE and TASTE were performed first with all participants and then again\nwith a reduced number of participants, excluding those who rated chocolate milk as less satisfying than water.\nExcluding participants did not alter the main connectivity results, though it slightly shifted obtained posterior\nprobability values. These values are presented below:\nModulatory influence of high-stress neutral and rewarding cues and taste on DCM connectivity parameters in\nthe full and reduced samples\nCovariates Nodes MAP Pp\nHigh-Stress Cue (n = 44) mean INS --> BNST -0.45 0.94\nmean BNST --> INS -0.34 0.91\nmean BNST --> OFC -0.45 0.99\nmean BNST --> NAcc -0.28 0.87\ndelta stress rating BNST --> OFC 0.19 0.83\nHigh-Stress Chocolate Milk Cues mean NAcc --> BNST 0.68 1\n(n = 44) mean BNST --> INS -0.28 0.86\ndelta stress rating n/a\nHigh-Stress Chocolate Milk Cuesa\nmean NAcc --> BNST 0.63 0.99\n(reduced participants, n = 38)\nmean BNST --> INS -0.21 0.68\ndelta stress rating n/a\nHigh-Stress Taste (n = 47) mean INS --> BNST -0.33 0.93\nmean NAcc --> BNST 0.28 0.65\ndelta stress rating BNST --> INS -0.11 0.58\nHigh-Stress Chocolate Milk Taste mean OFC --> BNST 0.18 0.54\n(n = 47) delta stress rating n/a n/a n/a\nHigh-Stress Chocolate Milk Tasteb mean n/a n/a n/a\n(reduced participants, n = 39) delta stress rating n/a n/a n/a\nBNST, bed nucleus of the stria terminalis; INS, insula; NAcc, nucleus accumbens; OFC, orbitofrontal cortex. MAP,\nmaximum a posteriori estimate; Pp, posterior probability. The Pp threshold cutoff was set at 0.75 (moderately\nstrong evidence). Continuous-time modulation is framed in terms of rates (measured in Hz), whereby a task (e.g.\nHigh Stress Cues) modulates the degree to which activity of an efferent node affects the continuous rate at which\nthe state of afferent node changes. aDCM for High-Stress Chocolate Milk Cue, excluding participants who rated\nchocolate milk as less satisfying than water (n = 38). bDCM for High-Stress Chocolate Milk Taste, excluding\nparticipants who rated chocolate milk as less satisfying than water (n = 39).\nAll participants provided responses for beverage, hunger and stress ratings across the beverage task stress\nconditions; thus, missing responses were not an issue. This is now specified in the methods section\nMS, p. 11:\n“All participants provided responses for beverage, hunger and stress ratings across both beverage task runs.”\n- For Figure 4 could the authors please also include the graph for the high stress water cue and high stress\nwater tasted so that the audience can see how those conditions differed from the milk conditions?\nResponse:\nWe thank the reviewer for the suggestion. Due to the high number of panels already present in Figure 4, we are\nconcerned that adding additional panels to depict the results for the high stress water cue and high stress water\nwould make the figure overcrowded and decrease readability of the main results of this work. However, we fully\nagree that it is of interest for readers to have the results of the high stress water cue and high stress water taste\nand these are now described in the MS and fully described in Table S4. Water taste did not result in any\nsignificant changes in effective connectivity during high stress.\nTable S4 – Modulatory influence of high-stress water cues and taste on DCM connectivity parameters\nNodes MAP Pp\nHigh Stress Water Cue INS --> BNST -0.44 0.81\nBNST --> OFC -0.37 0.89\nHigh Stress Water Taste - - -\nBNST, bed nucleus of the stria terminalis; MAP, maximum a posteriori estimate; INS, insula; OFC, orbitofrontal\ncortex. Pp, posterior probability. The Pp threshold cutoff was set at 0.75 (moderately strong evidence).\nContinuous-time modulation is framed in terms of rates (measured in Hz), whereby a task (e.g. High Stress\nWater Cues) modulates the degree to which activity of an efferent node affects the continuous rate at which\nthe state of afferent node changes.\n- The authors posit that inhibition of the NAcc via the BNST could promote reward-seeking behaviors to\n“ameliorate this striatal dampening” that happens during chronic stress. The reviewer thinks this\ninterpretation is over-stated. While the psychosocial stress task used in the study results in dampened\nNAcc activity, these data are from human neuroimaging studies, where we are unable to confirm the\nsource or underlying biological mechanism of the dampened activity. Can the authors provide a more\ndigestible interpretation of these data that align with the level of measurement we are able to achieve in\nhuman neuroimaging studies?\nResponse:\nWe thank the reviewer for raising this important point. We agree that our original interpretation could overstep\nthe bounds of inference permitted by human fMRI. In response, we have revised the relevant section of the\ndiscussion to better reflect the limitations inherent to fMRI-based effective connectivity modeling and to avoid\nimplying causality or directionality that cannot be empirically confirmed. We have rephrased the sentence in\nquestion and tempered the language used in the discussion section regarding our obtained findings.\nMS, p. 8:\n\"This pattern of BNST-to-NAcc inhibition may reflect stress-induced modulation of striatal responsivity, consistent\nwith preclinical models implicating this circuit in the regulation of reward-seeking behaviour under stress.\nHowever, given the limitations of human neuroimaging in measuring cellular-level mechanisms, this\ninterpretation is still speculative and grounded in broader preclinical literature on BNST-NAcc circuits.\"\nMinor concerns:\n- Please spell out HS_Cue, etc in figures and text to make it easier for readers to understand the conditions.\nResponse:\nWe thank the reviewer for the suggestion. We have replaced all instances of task-related acronyms with the full\nname in all tables, figures and main text, for easier readability.\n- The reviewer appreciates the authors performing LOO-CV and utilizing the data to also perform predict\nanalyses related to the relation between stress scores and brain effective connectivity. Could the authors\nreorganize the methods to reflect the organization of the Results section?\nResponse:\nWe have rearranged the methods section to reflect the organization of the results section as suggested.\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the\nNature Communications initiative to facilitate training in peer review and to provide appropriate recognition\nfor Early Career Researchers who co-review manuscripts.\nResponse:\nWe thank the Reviewer for the note. We appreciate the time and input that the Reviewer contributed to the\nreview process.\nIn this interesting and well-written report, the authors examined effective connectivity between the bed nucleus\nof the stria terminalis (BNST) and regions involved in reward, goal-directed decision-making, gustatory\nprocessing, and interoception during a beverage-tasting task implemented in the context of a brief stressor. Data\nobtained on a 7T scanner were analyzed using dynamic causal modeling and parametric empirical Bayes, with\nfindings reported as stress-induced inhibitory influence of BNST on cortico-striatal regions (nucleus accumbens,\norbitofrontal cortex, anterior insula) during cue exposure and beverage taste. There are several strengths of the\nwork, including a well-reasoned set of hypotheses, a focus on the BNST made possible through inclusion of\nhigh-resolution data from 7T, a relatively large sample size of healthy control subjects for this type of task\n(gustatory processing; n=~48), and rigorous implementation of DCM modeling and PEB approaches. That said,\nthe premise of the study rests on how the study was designed with respect to the stress induction, which\nunfortunately was methodologically weak. There are also challenges to understanding the specificity of the\nfindings in relation to the stress induction and to rewarding stimuli, due to the analytic approach. Together, these\nconcerns call into question the primary conclusions of the paper. Comments to detail these and related concerns\nare stated below.\nResponse:\nWe thank the Reviewer for their thoughtful and constructive summary, and for highlighting several key strengths\nof our study, including the rationale, use of high-resolution 7T imaging, and our relatively large sample.\nWe appreciate the Reviewer’s concerns regarding the stress induction and the specificity of our findings and have\ntaken care to address these in the revised manuscript by clarifying the experimental design, providing additional\njustification and context for the stress manipulation, and highlighting limitations wherever relevant.\n1) The study design included a standard gustatory (beverage) cue and taste task implemented in the scanner in\ntwo runs. Before the first run of the beverage task, a “low-stress” induction was completed (in the scanner),\ninvolving 2 minutes of counting backwards from a whole number. Before the second run, a “high-stress” induction\nwas completed (also in the scanner), involving 2 minutes of serial subtractions; errors or responding “too slowly”\nrequired starting over. There are several issues with this task and how the task was implemented.\na. The description of the task does not state any of the following critical details for the “high-stress” task:\n(a) from what number the subtraction started, (b) what numbers were subtracted serially, (c) who gave the\nfeedback, and how this was given (i.e., could the participant see the person giving feedback?), (d) what “too slow”\nmeant in terms of giving feedback on response speed. All of these details are important for interpreting results.\nResponse:\nWe thank the Reviewer for highlighting the need for greater detail in the description of the high-stress task.\nWe have revised the manuscript and Supporting Information to specify (a) that for the high-stress manipulation\nparticipants began subtracting from 1022, (b) that the step size was 13, (c) that feedback was delivered via\nintercom while participants viewed a two-minute timer on an MRI-compatible screen, and (d) that negative\nfeedback on response speed was delivered during brief pauses or delays in participant’s verbal responses. The\ntiming threshold for this feedback was based on the participant’s initial response tempo and applied consistently\nto simulate time pressure. All participants received negative feedback and were told that their performance was\n\"below average\" regardless of their actual performance. These clarifications are now reflected in the main text\nand described in further detail in the Supporting Information, and reproduced below:\nSI, p. 3–4:\nParticipants completed two structured stress induction tasks administered prior to each beverage fMRI task\nsession. The first induction (low-stress condition) involved an easy, two-minute serial subtraction task with no\n26\nevaluative feedback, while the second (high-stress condition) incorporated a difficult, two-minute serial\nsubtraction task with negative social evaluative feedback, modelled on the Montreal Imaging Stress Task (MIST).\nIn both conditions, participants viewed an on-screen digital clock (a simple rotating arm which remained white in\ncolour for the first minute and turned red for the second minute). In both conditions, participants reported\nresponses aloud and interacted with one researcher via an intercom. Below, to facilitate replication, we\nreproduce the relevant aspects of the protocol.\nLow-stress condition (Session 1)\n1. Researcher speaks via intercom: “Now we are going to start the beverage task. But before, we would like\nyou to count backwards from 1022 out loud. For example, 1022, 1021, 1020, et cetera. If you make a\nmistake, don’t worry, you can just start over from 1022 again. There will be a simple clock on the screen\nshowing passage of seconds, but take your time, we will tell you when to stop.\n2. The task begins when the participant is ready. After 30 seconds, researcher speaks (via intercom): “Now,\nplease count backwards from 100”.\n3. After 30 seconds, researcher speaks: “Thanks. Now, could you please count backwards from 50”.\n4. After 30 seconds, researcher speaks: “Thanks. Now, could you please count upwards, from 10 to 40 (for\nexample, 10,11,12,13, et cetera)”.\n5. If the participant completes this task with time remaining, researcher speaks: “Could you please count\nupwards from 70 to 100”.\nHigh-stress condition (Session 2)\nThis condition employed a modified version of the MIST. Participants were instructed to subtract 13 serially from\n1022 (and the resulting number), and to verbally report each result aloud under time pressure. Performance was\ncontinuously evaluated by the experimenter, with participants required to restart from 1022 upon errors or\ndelayed responses. Regardless of whether errors were made, participants received negative evaluative\nfeedback. The procedure is described below:\n1. Researcher speaks via intercom: “You will be asked to subtract the number 13—one three—from 1022,\nthen subtract 13 again from the new number—the result—and continue doing this out loud... So, 1022\nminus 13 and then subtract 13 from the result, and so on and so forth. We will be recording your answers\nand timing you. If you make a mistake or take too long to say the next number, you will be asked to start\nover from 1022. Like the previous task, there will be a simple clock on the screen. We will also be\ncomparing your response time to others.\n2. The task begins when the participant is ready. Participants receive immediate negative feedback during\nthe task if they make an error and are asked to restart—\"That was incorrect, please start again from\n1022”—and, if the participant is slow to respond (judged against the participant’s response tempo in the\ncontrol task), the researcher states: \"Your response is too slow, please start again from 1022”. Here, the\nimportant thing is that the researcher exploits pauses to insert negative feedback.\n3. At the end of the two minutes, researcher states: \"Thank you. Your response time was below average\".\n27\nb. Stressors that are 2 minutes in length are not convincing in terms of actually inducing stress, and the task\ndoes not appear to include components of uncontrollability/unpredictability. The negative evaluative threat\ncomponent is not particularly convincing either. Added to this, objective measures of stress either were not\ncollected (i.e., serum or salivary cortisol) or in the case of electrodermal activity (EDA) did not show any\ndifference between “low-stress” and “high-stress” (p = 0.74). A large body of research on robust psychosocial\nstress tasks, such as the Trier Social Stress Test, the Maastricht Acute Stress Task, and the Montreal Imaging\nStress Task are much longer (15-20 minutes), and studies which employ them almost always include an objective\nmeasure of HPA-axis reactivity and/or include a well-matched control version which is not subject to order effects\n(as was the case here, where all subjects first had “low-stress” then “high-stress” stress conditions).\nResponse:\nWe appreciate the Reviewer’s concerns regarding the structure of our stress manipulation. We offer the following\npoints of clarification and justification:\nDuration of stressor\nWe note that brief in-scanner manipulations of similar duration have been used in prior neuroimaging studies\n(for review, see Berretz et al., 2021). Indeed, the Montreal Imaging Stress Task (MIST) — which the Reviewer\ncites as an example of a robust psychosocial stress paradigm — can involve stress blocks as short as 2 minutes,\ndepending on the participant’s performance during a preceding calibration session (see Dedovic, 2005). Our\nown design followed these precedents, balancing the need for stress induction with the practical constraints of\na neuroimaging environment, particularly in the context of ultra-high-field (7T) MRI. We now acknowledge the\nduration of stressor as a potential limitation in the Discussion:\nMS, p. 8:\n“Finally, the stress induction period was relatively brief (two minutes), and although we observed both subjective\nand physiological signs of stress, neuroendocrine measures were not collected to provide additional confirmation\nof hypothalamic-pituitary-adrenal-axis engagement.”\nSee:\n• Berretz, Gesa, et al. \"The brain under stress—A systematic review and activation likelihood estimation\nmeta-analysis of changes in BOLD signal associated with acute stress exposure.\" Neuroscience &\nBiobehavioral Reviews 124 (2021): 89-99.\n• Dedovic, Katarina, et al. \"The Montreal Imaging Stress Task: using functional imaging to investigate the\neffects of perceiving and processing psychosocial stress in the human brain.\" Journal of Psychiatry and\nNeuroscience 30.5 (2005): 319-325.\nUnpredictability and uncontrollability\nThe high-stress task incorporated pseudo-contingent negative feedback based on participants’ response tempo,\nrather than accuracy alone. This design introduced an element of unpredictability, as participants could not\nreliably determine the performance criteria and were consistently told to restart based on their accuracy and\nspeed responding. This created a sense of negative social evaluation and uncontrollability, all hallmarks of\npsychosocial stress (Kemeny, 2023). This was behaviourally endorsed by participants in the ratings they\nprovided after the low-stress and high-stress induction (paired t-test: t = 12.41, p < .0001, n = 48, η²p = 0.7662;\n(47)\nmean difference between low and high-stress induction was 2.771 [SD of differences = 1.547, 95% CI = 2.322-\n3.220, Cohen’s d = 1.7912]). We have clarified this in both the main text and Supporting Information.\nSee:\n28\n• Kemeny, M. E. The Psychobiology of Stress. Curr Dir Psychol Sci 12, 124–129 (2003).\nPhysiological evidence of stress\nAlthough electrodermal activity (EDA) did not significantly differ between conditions, we now show that the high-\nstress condition was associated with reductions in HRV, consistent with sympathetic arousal. As expected, the\naverage root mean square of the successive differences (RMSSD) in pulse signals—a measure of HRV—during\nthe start of the high-stress run was lower than during the low-stress run, t(43) = -2.25, p = 0.029.\nWe also note that stress, as inferred via HRV, was sustained after the high-stress task (relative to the low-stress\ntask). Specifically, we evaluated the sustained after-effects of high stress induction by analysing the complete\ntime course of RMSSD using a sliding-window approach (10-second windows stepped every 2 seconds) across\nthe beverage task. We then segmented the data into early, middle, and late task periods (first, middle, and final\n20 seconds) to evaluate trends across time. The within-participants repeated-measures ANOVA showed a\nsignificant Condition × Time interaction F(2, 86) = 3.20, p = 0.045, indicating that RMSSD evolved differently\nover time between the low- and high-stress conditions. Additional paired t-test comparisons showed that RMSSD\nincreased significantly more from early to middle in the LS condition compared to the high-stress condition, t(43)\n= 2.25, p = 0.029. Similarly, from early to late, RMSSD increased in the low-stress condition but decreased in\nthe high-stress condition, t(43) = 2.03, p = 0.048. These results are now all presented in the manuscript.\nNo cortisol sampling\nWe agree that salivary cortisol and other HPA-axis-related markers would offer additional confirmation of stress\ninduction, but these were not feasible in the current design. We have now noted this as a limitation (see\nreproduced text above).\nOrder effects\nWe employed a fixed order of stress conditions (low-stress followed by high-stress) to minimize carryover effects\n(prolonged physiological arousal, for example). Our design thus prioritized preserving an unbiased baseline for\nthe low-stress condition. This approach is consistent with precedent in neuroimaging protocols; however, we\nacknowledge that counterbalancing stress order could be implemented and recommend it as a future design consideration.\nSee:\n• Wang, J., et al., 2005. Perfusion functional MRI reveals cerebral blood flow pattern under psychological\nstress. Proceedings of the National Academy of Sciences.\n• Wheelock, et al., 2016. Prefrontal cortex activity is associated with biobehavioral components of the\nstress response. Front. Hum. Neurosci.\nMS, p. 8:\n“Furthermore, we cannot exclude potential order effects with our study design (low-stress always preceded the\nhigh-stress condition), which we implemented to minimize carry-over effects as in previous work. Future\nstudies should aim to implement washout periods or counterbalanced designs to better isolate stress-\nspecific contributions.\nc. Given the design that did not adequately address the control condition issue, it is nearly impossible to ascertain\nwhether results stemming from the “high-stress” condition were simply a reflection of being in the scanner for as\nlong as they were at that point, having to complete a somewhat challenging task involving cues, water/chocolate\nmilk tastes, swallowing, rinsing, etc.\nWe acknowledge that we cannot fully attribute effects to the stress induction alone. However, the high-stress\ncondition was associated with sustained reductions in HRV and elevated self-reported stress relative to the low-\nstress condition, despite both involving identical sensory and motor demands. These findings suggest that the\nobserved effects are unlikely to reflect general task fatigue or time-in-scanner effects.\nd. The references provided as suggesting this task has been validated are insufficiently rigorous as examples\nfor the current study. The first reference (77; Kistenmacher et al) used a very different cross-over design, in\nwhich subjects completed 2 visits – one that was control/no stress and one that was stress, and the stress\ninduction was implemented for a full 5 minutes. In addition to these differences, cortisol and ACTH were collected\nand analyzed to demonstrate a clear stress response in participants during the stress induction visit. The task\nused in the current study does not match this design.\nResponse:\nWe thank the Reviewer for their close reading. We agree that the design reported in Kistenmacher et al.’s\npublication differs in structure (e.g., crossover, inclusion of endocrine measures), and we now clarify in the\nmanuscript that our reference to prior work was intended to establish precedent for using brief, in-scanner math-\nrelated stressors, not to imply use of an identical methodology. With regard to clear indications of stress,\nwe concede that the clearest evidence for stress induction in our study comes from self-reported ratings\n(primarily) and autonomic HRV markers (secondarily):\nMS, p. 10:\n“There is a strong precedent for using brief, in-scanner math-related tasks to induce stress in the context of\nfMRI76–79. “\ne. The other references listed do not provide sufficient detail in their methods to match the current study either,\nand also did not include objective measures of HPA-axis reactivity, so it is challenging to accept these as\n“validation” of the task used in the current study.\nResponse:\nWe acknowledge the Reviewer’s concern regarding the absence of HPA-axis measures in the cited studies.\nHowever, we emphasize that our use of these references was not to claim identical validation, but to highlight\nconverging support for brief verbal tasks that elicit stress in MRI environments. While endocrine markers can be\nvaluable, we view them as complementary rather than essential for establishing task validity in a neuroimaging\ncontext.\nOur approach demonstrates face validity, statistical conclusion validity (with significant effects in the\nhypothesized directions), and convergent validity with prior findings using related paradigms. Moreover, the\nobserved changes in self-reported stress ratings and HRV provide multiple lines of evidence for successful stress\ninduction.\nf. Although the report does provide analysis of self-reported “stress” ratings, again, given the design, it is not\npossible to know whether the differences in ratings reported after the “low-stress” induction compared to those\nreported after the “high-stress” induction simply reflect the influence of demand characteristics, because there is\nno control arm/visit. This is particularly likely to be the case given the absence of objective, physiological\nmeasures showing a stress response (i.e., cortisol/ACTH were not measured, and EDA response showed no\ndifference between conditions). Unfortunately, based on the above concerns regarding the design and\nimplementation of the stress induction, the discussion and conclusions are not well supported by the data.\nResponse:\nWe appreciate the Reviewer’s concern regarding the interpretability of self-reported stress ratings in the absence\nof a control arm or endocrine markers. However, we respectfully challenge the assumption that the absence of\ncortisol or ACTH measurements precludes meaningful inferences about experienced stress. Psychological\nstress is fundamentally a subjective, appraisal-based state, and as such, self-reports remain the most direct and\ninterpretable measure of that construct (consider that, in a related area of research, there is no evidence for\nunaware fear conditioning).\nSee:\n• Boddez, Y, et al. \"Rating data are underrated: Validity of US expectancy in human fear conditioning.\"\nJournal of behavior therapy and experimental psychiatry 44.2 (2013): 201-206.\nWhile physiological markers can offer converging evidence, they are not specific to stress, nor are they\nuniversally sensitive to stress. Indeed, research in threat processing shows that physiological responses (such\nas skin conductance, heart rate or fear potentiated startle) are contingent on subjective awareness of the\nstressor—there is no reliable evidence for contingency-unaware stress or fear responses at the physiological\nlevel. This further supports our focus on stress ratings.\nSee:\n• Mertens, G, and Iris M. Engelhard. \"A systematic review and meta-analysis of the evidence for unaware\nfear conditioning.\" Neuroscience & Biobehavioral Reviews 108 (2020): 254-268.\nWe also reiterate that our stress induction produced consistent differences across multiple response modalities,\nincluding self-report ratings and HRV, supporting the validity of the manipulation despite the absence of HPA-\naxis sampling.\n31\n2) The authors mention that EDA values were “higher after high-stress induction compared to the low-stress\ninduction, though this difference did not reach significance”. This is not only evidence that the stress induction\ndid not actually invoke a physiological response, but is also misleading to state. The mean difference stated was\n-0.0006 (p = 0.74); this is negligible and should not be framed as “higher after high-stress”, even if it didn’t reach\nsignificance. Instead, the report of the finding should be revised to make it clear that EDA was not different\nbetween conditions, without further elaboration/interpretation.\nResponse:\nWe thank the Reviewer for raising this point. As noted in our response to Reviewer 1, we have revised the text\nto clearly state that there was no significant difference in EDA between conditions, and we have removed any\nlanguage implying a trend or directional effect.\nMS, p. 5:\nEDA values were not significantly higher after high-stress induction compared to the low-stress induction (t = -\n27\n0.31, p = .76, n = 28, mean difference = -0.006 [-0.05, 0.04]; see Figure S1).\nSI, Fig. S1:\nDifferences between conditions did not reach statistical significance (t = -0.31, p = .76, n = 28, mean difference\n27\n= -0.006 [-0.05, 0.04]).\n3) As it is described, it appears that the DCM model examines beverage cues (chocolate milk + water, or\nchocolate milk alone) and beverage tastes (again, chocolate milk + water, or chocolate milk alone) specifically\nin the “high-stress” condition, without comparing to the “low-stress” condition. Therefore, any findings related to\nthe analytic approach prevent inferences related to specificity to stress (especially given concerns described\nunder point 1) and could alternatively have been interpreted as findings specific to completing run 2 of the\nbeverage taste task, not to stress.\nResponse:\nOur analysis focuses on characterizing effective connectivity during the high-stress condition as compared to\nthe low-stress condition. Our DCM specification included concatenated time series across both task runs to allow\nfor this analysis. This design allowed for the investigation of network dynamics specifically within cue\npresentation and taste processing during a high-stress state compared to the low-stress state. We have provided\nconverging evidence—including reduced HRV and elevated self-reported stress—that the stress induction was\nsuccessful. We have clarified this aim in the revised manuscript to avoid any confusion regarding the intended\nscope of inference.\nMS, p. 7:\nHere, we characterize effective connectivity during acute stress and show that exposure to beverage cues and\ntaste receipt is associated with BNST-driven inhibition of cortico-striatal regions.\n4) Similarly, the model does not appear to contrast chocolate milk to water, in some cases these are combined\n(models 1 and 3), reducing confidence that the rewarding aspect of the tastants is involved, yet much of the\ndiscussion of these findings is framed around reward. This is confusing and perhaps misleading given lack of\nspecificity to the rewarding stimulus (chocolate milk), which was rated as significantly more pleasant than the\nwater stimulus.\nResponse:\nIn contrast to standard GLM analyses, DCM does not support direct subtraction-based contrasts (e.g., chocolate\nmilk vs. water), as it models condition-dependent changes in connectivity rather than activation differences. Our\nprimary aim was to examine effective connectivity modulatory effects during a high-stress state involving cues\nand tastes relative to a low-stress stress. To address the issue of reward specificity and in agreement with the\nother Reviewer’s suggestion, we conducted a separate analysis using water-only trials and present the key\nresults below.\nTable S4. Modulatory influence of high-stress water cues and taste on DCM connectivity parameters\nNodes MAP Pp\nHigh Stress Water Cue INS --> BNST -0.44 0.81\nBNST --> OFC -0.37 0.89\nHigh Stress Water Taste - - -\nBNST, bed nucleus of the stria terminalis; MAP, maximum a posteriori estimate; INS, insula; OFC, orbitofrontal\ncortex. Pp, posterior probability. The Pp threshold cutoff was set at 0.75 (moderately strong evidence).\nContinuous-time modulation is framed in terms of rates (measured in Hz), whereby a task (e.g. High Stress\nWater Cues) modulates the degree to which activity of an efferent node affects the continuous rate at which\nthe state of afferent node changes.\n5) A minor point to note is that the authors report excluding 6 subjects whose chocolate milk pleasantness ratings\nwere lower than the water pleasantness ratings, presumably because their ratings would reflect that the\nchocolate milk was not pleasant, but they do not mention whether data were excluded for any subjects who rated\nstress as similar across runs. This should be considered because those who rated stress as similar across runs\n(i.e., those for whom stress ratings did not change) or were lower for the “high-stress” relative to “low-stress”\nwould similarly not fit their assumption that only data for those who found it to be stressful should be included.\nResponse:\nWe thank the Reviewer for their observation. We agree that intra-participant stress ratings across conditions is\nan important consideration. No participant rated the low stress induction as more stressful than the high stress\ncondition. A paired-samples t-test confirmed that perceived stress was significantly higher in the high-stress\ncondition (t(47) = 12.41, p < .0001, n = 48, η²ₚ = 0.7662).\nThe mean difference in stress ratings was 2.77 (SD =\n1.55), with a 95% confidence interval of [2.32, 3.22], and a large effect size (Cohen’s d = 1.79), indicating a\nrobust manipulation.\nTitle:\nBed nucleus of the stria terminalis connectivity during food cue and taste processing under stress\nManuscript: NCOMMS-25-24841A\nFor convenience, all comments are provided in boxed format with responses below. Quoted text from the\nmanuscript and supporting information appears in blue. A clickable Table of Contents is provided below to\nfacilitate navigation.\nContents:\nReviewer 1 ......................................................................................................................................................... 2\nReviewer 2 ......................................................................................................................................................... 4\nReviewer 3 ....................................................................................................................................................... 12\nReviewer 4 ....................................................................................................................................................... 13\nReviewer 1\nI want to commend the authors for their largely convincing and thoughtful revision of the previously raised\nmajor issues. I also briefly reviewed the documentation of the analysis code and the associated material, and\nfound this to be very useful to establish the robustness of the revised analysis.\nAlthough I do not want to unnecessarily drag on the review process, I have one remaining major concern that\nI believe is important. The revised manuscript now presents one HRV measure to \"confirm\" the differences\nbetween the stress conditions. However, there is no good rationale to use this single measure, specifically\nsince many other stress fMRI studies use heart/pulse rate, and not only one measure of HRV (which are\nthought to reflect partly separate processes of autonomic balance during stress exposure).\nIn my opinion, it should be transparently reported that EDA was the primary objective stress measure and\nthat those HR(V) analyses are post hoc now, and the revised report should include HR and, ideally, additional\nHRV measures. These indices do not necessarily have to be interpreted in isolation, but the question is now\nwhether RMSSD was cherry-picked or not (there is considerable evidence concerning such publication\nbiases in the HRV literature, so the paper should not follow these examples).\nResponse:\nWe thank the Reviewer for commending our revision and for raising these important and points. We\nacknowledge that using the term “confirm” was imprudent and this has now been altered in the manuscript.\nIn our revised manuscript, we now indicate that analyses of heart rate variability (HRV) were post hoc. As\nsuggested, we have re-examined HRV across task runs using multiple complementary measures.\nSpecifically, tests for short-term variability measures (pNN20 and pNN50) and the nonlinear Poincaré ratio\n(SD1/SD2), assessing the relationship between short- and long-term heart-rate fluctuations, were run.\nHigh-stress runs showed lower pNN20 (t = 2.96, p = 0.0048) and pNN50 (t = 2.39, p = 0.021), a higher\n(47) (47)\nlow- to high-frequency ratio (t = −2.61, p = 0.012), and a reduced SD1/SD2 ratio (t = 2.07, p = 0.044).\n(47) (47)\nThese converging results indicate reduced parasympathetic modulation and a relative increase in\nsympathetic activation during the high stress run, consistent with greater physiological arousal. These\nfindings are now reported in the results section and additional information on these analyses can be found in\nthe methods section, supplementary information and on the manuscript’s corresponding github page.\nMS, p. 5:\n“These results suggest that (1) the stress manipulation was successful in increasing participants’ self-\nreported stress levels […].”\nMS, p. 5:\n“Physiological stress was primarily indexed using tonic electrodermal activity (EDA; detailed in Methods),\ncomplemented by exploratory analyses and heart rate variability (HRV) indices.”\nMS, p.17:\n“In addition, successive inter-beat intervals were converted into difference vectors, from which the proportion\nof successive differences exceeding 20 ms (pNN20) and 50 ms (pNN50) was calculated. The nonlinear\nPoincaré index (SD1/SD2) was derived from the standard deviations of instantaneous beat-to-beat changes\n(SD1) and longer-term variability (SD2). Mean values from low- and high-stress runs were compared using\npaired-sample t-tests.”\nRemarks on code availability:\nI briefly reviewed the code and documentation to verify consistency with the report in the manuscript. I did\nnot try to rerun the analysis (and I truly believe reviewers need to be given more time to carefully review a\nrevision if this is intended)\nResponse:\nWe thank the Reviewer for reviewing the code and documentation. We appreciate acknowledgment of the\ncode’s consistency with the report.\nReviewer 2\nThe authors made substantial edits and added supplemental materials to the present manuscript\ninvestigating BNST effect connectivity during taste under stress. The clarification and revision of the regions\nof interest was appropriate. The addition of the analysis for the low stress and water conditions is helpful.\nHowever, that raises a new question, also raised by Reviewer 4 about including the control conditions in the\nmain analyses. Since the goal of the study is to examine human food cue and taste processing under stress,\nit is critical to isolate the effect of stress by including the control condition. In the rebuttal letter the authors\nstate that DCM does not support direct subtraction-based contrasts. However, there are several ways to\ninclude the ‘control’ conditions in the analysis. DCM does support having a modulatory input that represents\neach of the conditions in one model. Another alternative would be to include low stress water as the implicit\nbaseline in the first-level models. At the very least the authors could include the control models in the figures\nand main results so that the reader can at least infer the effective connectivity changes that are modulated\nby the tasks of interest.\nResponse:\nWe thank the Reviewer for raising this question regarding the inclusion of all four task conditions. In response,\nwe have now inverted all cue- and taste-related conditions in the subject dynamic causal models (DCMs). In\nthese models, the influence of all four conditions—low-stress water, low-stress milk, high-stress water, and\nhigh-stress milk—was modelled simultaneously. These revised models were specified such that all relevant\ncontrasts could be recovered.\nWith the A- and B-matrices specified exactly as in the main analyses, the expanded four-condition model did\nnot introduce new parameters. However, it did require the same set of B-parameters to explain a substantially\ngreater amount of condition-specific structure. As expected, this increased modelling demand attenuated the\nstimulus-related effects at the group level (under parametric empirical Bayes modelling and Bayesian model\naveraging), with the relevant parameters shrinking.\nImportantly, the modulatory influence of stress remained consistent across all model specifications.\nAs suggested by the reviewer, we have included the expanded four-condition models for the reader’s\nreference in the manuscript and report the complete results in an expanded table in the supporting\ninformation. We have retained the original models in the main text as the principal analysis as this is\nconsistent to our original study aims around specifically examining influence of stress on BNST connectivity.\nMS, p. 7:\n“Results from additional models—specified identically to the main analyses except that the modulatory inputs\nwere coded to represent all four conditions simultaneously—are also reported in the Supporting Information\n(Table S5).”\nSI, p. 13:\nModulation\nModel Connection LS-Water LS-Milk HS-Water HS-Milk\nCue BNST → NAcc 0.209 0.209 –0.209 –0.209\nBNST → OFC 0.321 0.321 –0.321 –0.321\nBNST → dmINS 0.226 0.226 –0.226 –0.226\nNAcc → BNST 0.228 0.228 –0.228 –0.228\nOFC → BNST 0.120 0.120 –0.120 –0.120\ndmINS → BNST 0.197 0.197 –0.197 –0.197\nTaste BNST → NAcc 0.113 0.113 –0.113 –0.113\nBNST → dmINS 0.193 0.193 –0.193 –0.193\nOFC → BNST 0.259 0.259 –0.259 –0.259\ndmINS → BNST 0.186 0.186 –0.186 –0.186\nThis table reports condition-specific modulation of effective connectivity (Hz) for each connection that\nexhibited a supra-threshold modulation at the group level (posterior probability Pp ≥ 0.75). All four\nconditions—Low-Stress Water, Low-Stress Chocolate Milk, High-Stress Water, and High-Stress\nChocolate Milk—were modelled simultaneously using two orthogonal modulators (Stress: +0.5/–0.5;\nStimulus: +0.5/–0.5), and the values reported here reflect the linear combination of estimates for each\ncondition. As the Stimulus modulator did not survive the Pp ≥ 0.75 criterion in either the cue or taste\nmodel, its coefficients were effectively zero at the group level, resulting in identical values for LS-Water\nand LS-Milk (and likewise for HS-Water and HS-Milk). Region abbreviations: BNST, bed nucleus of the\nstria terminalis; NAcc, nucleus accumbens; OFC, orbitofrontal cortex; dmINS, dorsal mid-insular cortex.\nOther major concerns:\nUse of causative language is still present throughout the manuscript. For example, the abstract, hypothesis,\nand DCM methods section use describes the connectivity changes as “upregulates” or “downregulates”. It\nwould be more accurate to say decreases/increases the strength or weakens/strengthens the effective\nconnectivity. DCM is a unique method that allows us to gather directional evidence in humans, but cannot\nprovide any evidence about up or down regulation.\nResponse:\nWe acknowledge the Reviewers’ concern and agree that there is a risk that these terms — “upregulation”\nand “downregulation” — evoke biological processes beyond the scope of the current models and analyses.\nThat said, with reference to the models’ modulatory parameters, these terms are suggested by the developers\nof DCM (Zeidman et al., 2019), and causal inferences are a key aspect of the causal (control-theoretic)\nperspective that motivates these models (Valdes-Sosa et al., 2011). In response, we have modified text at\nthe end of the Introduction to clarify the use of this terminology in relation to the (mean-field) activity of\nneuronal populations. Below, we reproduce some of the key adjustments.\nReferences:\nZeidman, P. et al. A guide to group effective connectivity analysis, part 1: First level analysis with DCM for\nfMRI. NeuroImage 200, 174–190 (2019).\nValdes-Sosa, P. A., Roebroeck, A., Daunizeau, J. & Friston, K. Effective connectivity: Influence, causality\nand biophysical modeling. NeuroImage 58, 339–361 (2011).\nMS, p. 3:\n“However, understanding food cue-related interactions within brain networks has been hampered by the\nfield’s reliance on correlational measures (i.e., functional connectivity), rather than methods that can infer\neffective connectivity—the signed, directed influence that one brain region exerts over another43, where a\npositive (negative) connection reflects a net excitatory (inhibitory) influence at the level of neuronal\npopulations. Though not superior to simpler functional connectivity models, dynamic causal modelling (DCM)\noffers a biophysically grounded framework for inferring effective connectivity, and has proven valuable in\naddressing such questions 44.\nHere, we leveraged ultra-high field 7-Tesla fMRI, DCM and controlled taste delivery via a gustometer to map\nthe effective connectivity between the BNST and key regions implicated in food cue and taste processing:\nthe NAcc, insula, and OFC. We hypothesized that food cue and taste processing under high stress would\ndownregulate (decrease the values of) BNST connectivity with other network nodes.”\nIn the revised discussion, the authors cite Flook et al for evidence of structural connectivity between the\nBNST and the mid-insula subdivision. In that paper, it appears that the only significant connectivity is with the\nanterior-ventral region of the insula, with very low connectivity with the mid-insula subregions. It might\nstrengthen the authors’ argument to provide support for BNST-insula connectivity using other methods, like\nresting state connectivity.\nResponse:\nWe thank the reviewer for this helpful suggestion. We now reference “Resting state connectivity of the bed\nnucleus of the stria terminalis at ultra-high field” (Torrisi et al., 2015), which provides complementary evidence\nfor BNST–insula connectivity using resting-state fMRI. This strengthens the anatomical and functional basis\nfor including this pathway in our model.\nMS, p. 3:\n“Evidence from diffusion-weighted MRI, resting-state fMRI and post-mortem studies have confirmed\nbidirectional white matter tracts between the BNST, NAcc, insula and OFC in humans28–32, such that these\nregions form an anatomical network.”\nReferences:\nTorrisi, S. et al. Resting state connectivity of the bed nucleus of the stria terminalis at ultra‐high field. Human\nBrain Mapping 36, 4076–4088 (2015).\nConsidering the authors found an effect of BNST-OFC effective connectivity on changes in stress scores, it\nwould be helpful to include a scatterplot showing the correlation between the two measures, in addition to\nthe predicted and actual scores graphs.\nResponse:\nThanks for the suggestion. We have now added a scatterplot illustrating the relationship between the subject-\nlevel connectivity parameters and stress sensitivity scores. This is now depicted in Figure S5.\nMS, p. 16:\n“An additional subject-level analysis (Figure S5) showed that individual differences in stress-evoked BNST-\nto-OFC modulation (encoded in DCMs) was positively associated with changes in subjective stress, providing\na complementary, non-PEB perspective on the same effect.”\nSI, p. 23:\nEach point reflects a single subject’s first-level estimate of high-stress-related modulation on the BNST-to-\nOFC connection (i.e., the individual B-parameter prior to entering the group-level PEB model), plotted against\nthat subject’s within-person change in subjective stress (z-scored). Because these are first-level estimates,\nthe values are statistically independent (no empirical Bayes shrinkage has been applied). The fitted line\nshows the ordinary least-squares regression line corresponding to the observed Pearson correlation (r =\n0.30, p = 0.44), illustrating a positive association between change in subjective stress and BNST-to-OFC\nmodulation.\nMinor concerns:\nThe authors clarified in the rebuttal that the cluster threshold (k = 10) was only used for visualization purposes.\nBut the methods in the manuscript still says, “A false discovery rate (FDR) corrected statistical threshold was\napplied (PFDR<.05) to all whole-brain GLM analyses with a 10-voxel cluster-extent threshold (KE ≥10\nvoxels)”.\nResponse:\nWe are happy to clarify this point. FDR peakwise set a voxelwise t-value cutoff (pFDR<.05) for each contrast,\nand this was our inferential threshold for the GLM analyses. The k≥10 criterion was added as a reporting\nrule, and we did not depict in the figures or describe in the tables any FDR-significant voxels that fell in\nclusters smaller than 10 voxels. This approach is consistent with methodological work showing that voxelwise\nFDR/FWE control is statistically robust, whereas clusterwise-alone inference is more prone to inflated false\npositives (particularly for small clusters). As such, a minimum cluster size was applied only for reporting\nclarity, consistent with methodological recommendations to avoid overinterpreting isolated voxels. (Woo et\nal., 2014, NeuroImage; Eklund et al., 2016, PNAS). We have revised the Methods to state explicitly that\nsignificance was determined by voxelwise FDR, with only clusters ≥10 voxels visualized and reported. We\nnow specifically state in the methods section that this value is applied for both visualization and reporting\npurposes:\nMS, p. 13:\n“A KE ≥10 extent threshold was used for visualization and reporting purposes to avoid over-interpreting\nisolated voxels.”\nThe authors acknowledge that the area of activation is not anterior insula and have updated the term to\ninsula. However, the insula is large and heterogenous, with well-described subregions that have distinct\ncytoarchitectures, structural and functional connectivity, and functions. For example,\ndoi:10.1093/cercor/bhs065 uses functional parcellation to identify large sub-regions and includes a meta-\nanalytic coactivation study to infer functions. The authors do refer to the region as “mid-insula” in the\ndiscussion, but it would be helpful to provide justification for selecting the appropriate sub-region and use\nthat name throughout.\nResponse:\nWe thank the reviewer for bringing up this point and agree that, given the heterogeneity of the insula, our\ninclusion of our mid-insula peak warrants further justification. As suggested, we have changed the name of\nour insula peak to “dorsal mid-insula (dmIns)” through the manuscript.\nThis dmIns peak was included in DCM model space as it delivered the highest peak activation values within\nthe insula for both our cue and taste GLM analyses. Moreover, the dmIns has been found in recent 7T work\nto serve as a core node within the human interoceptive network. This complements other research by\nLuchsinger et al. demonstrating that the mid-insula-to-BNST pathway that is recruited during stress and is\nmodulated during subsequent anxiety-like behaviour. We now cover this distinction in the manuscript with\nthe following text:\nMS, p.9:\n“In response to an acute stressor, activity of dmINS-to-dorsal BNST corticotropin releasing factor (CRF)-\npositive neurons is associated with active struggling behaviour - a physical stress-coping mechanism in\nrodents64. Relatedly, dmINS has been found in recent 7T research to function as a primary interoceptive hub\ntightly coupled with subcortical allostatic centers regions, positioning it as a key node linking bodily state\nencoding with stress-responsive control.”\nZhang, J., Chen, D., Deming, P. et al. Cortical and subcortical mapping of the human allostatic–interoceptive\nsystem using 7 Tesla fMRI. Nat Neurosci 28, 2380–2391 (2025).\nLuchsinger, J.R., Fetterly, T.L., Williford, K.M. et al. Delineation of an insula-BNST circuit engaged by\nstruggling behavior that regulates avoidance in mice. Nat Commun 12, 3561 (2021).\nClarify in the methods that the DSM-IV/ICD-10 version of the MINI was used for screening.\nResponse:\nWe have made this change to the manuscript\nWhat was the number of left vs. right strongest task-evoked peaks for the OFC and insula? It looks like all\npeaks were from one side (right?). If this is the case, it would be helpful to include a discussion of laterality.\nResponse:\nThe reviewer is correct that all the strongest peaks were obtained from the right hemisphere for the OFC and\nthe insula, though significant activation peaks were also detected for the cue and taste conditions in the left\nhemisphere. We have clarified this point in the methods section with the following sentence:\nMS, p. 14:\n\"The OFC and dmINS were modelled unilaterally, using the strongest task-evoked peaks from the GLM, in\nline with DCM guidelines prioritizing dominant hemispheric activation, which in this sample were located in\nthe right hemisphere.”\nNotably, both the OFC and dmINS nodes of our model showed right-lateralised task-evoked peaks,\nconsistent with work indicating preferential engagement of right insula in interoceptive/salience and stress-\nrelated processing and asymmetric reward-related signalling in right OFC.\nAs suggested, we have also included the following text around laterality in the discussion:\nMS, p. 8:\n“Moreover, our right-lateralised activation pattern aligns with evidence that interoceptive–salience\ncomputations and reward-valuation signals show stronger right-hemisphere engagement, indicating that our\nDCM nodes captured this broader functional asymmetry”\nLopez-Persem, Alizée et al. “Differential functional connectivity underlying asymmetric reward-related activity\nin human and nonhuman primates.” Proceedings of the National Academy of Sciences of the United States\nof America vol. 117,45 (2020): 28452-28462. doi:10.1073/pnas.2000759117\nAdamic, Emily M et al. “Hemispheric divergence of interoceptive processing across psychiatric disorders.”\neLife vol. 13 RP92820. 13 Nov. 2024, doi:10.7554/eLife.92820\nFigure S2 should have a legend in the figure for what the activation colors indicate.\nResponse:\nWe have now updated this figure.\nSI, p. 19:\nOverlaid statistical maps showing thresholded group-level activation (FDR-corrected, p < 0.05) for high-stress\ncue (shown in orange) and high-stress taste (shown in blue) conditions across the bed nucleus of stria\nterminalis (BNST), dorsal mid-insula (dmINS), nucleus accumbens (NAcc) and orbitofrontal cortex (OFC).\nROIs are shown on a standard brain surface in MNI space. Colour overlays indicate spatial extent and overlap\nof activation between tasks within each ROI. BNST Taste activation outlined in white, NAcc Taste activation\noutlined in red.\nThe authors should consider simplifying language, e.g. “By contrast, focusing on taste receipt under high\nstress, we found that taste receipt downregulated effective connectivity from the insula to BNST (Figure 4c);\nhowever, this effect did not exceed our imposed posterior probability threshold when considering chocolate\nmilk and water taste receipt separately.” This could be more succinctly stated as there was a main effect of\nhigh stress, regardless of beverage.\nResponse:\nWe thank the Reviewer for this suggestion. We have made a revision to the text to more clearly (and\nconcisely) describe the posterior inference.\nMS, p. 7:\n“By contrast, focusing on taste receipt under stress, we found that taste receipt downregulated insula-to-\nBNST effective connectivity (Figure 4c); however, this effect did not exceed our imposed posterior probability\nthreshold when considering chocolate milk and water taste receipt separately.”\nSimilarly, the authors had an additional analysis “treating stress sensitivity as a covariate”—can the authors\nbe clearer about this? If the covariate was difference between high and low stress induction, using “change\nin subjective stress” would be more precise in describing the strength of the stress induction for each\nparticipant. Stress sensitivity typically refers to a trait, which was not measured in this study. Also please\nclarify if the calculated score is based on average ratings for low and high stress or computed during specific\nperiods, e.g., late period of low-stress to early period of high stress.\nResponse:\nWe thank the Reviewer for flagging this lack of precision. All instances of “stress sensitivity” have been\nadjusted to the more descriptive “change in subjective stress”. Below, we reproduce changes made to the\nMethods section, which we hope makes the calculation of this variable clear.\nMS, p. 15:\n“One PEB model (Figure 4d) included the effect of change in subjective stress, quantified as the change\n(delta) from low- to high-stress induction46 (i.e., the self-reported stress score after low-stress induction minus\nthe self-reported stress score after high-stress induction).”\nReviewer 3\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the\nNature Communications initiative to facilitate training in peer review and to provide appropriate recognition\nfor Early Career Researchers who co-review manuscripts.\nResponse:\nWe thank the Reviewer for the note. We appreciate the time and input that the Reviewer contributed to the\nreview process.\nThe authors have responded thoughtfully to prior critiques, including providing details regarding the task\ninstructions for the stressor, adding HRV data as a physiological measure of stress, and noting in the\ndiscussion several limitations related to the design of the stress task. However, there remain unresolved\nconcerns, some of which require clarification and revisions to wording to soften interpretation, others which\nrequire significantly more detail and justification for analytic approach and justification for interpretation (HRV\ndata). Additionally, adherence to the SAGER guidelines appears inadequate and requires additional analysis\nand/or presentation of data disaggregated by sex.\nResponse:\nWe thank the Reviewer for acknowledging the revisions made in response to earlier feedback, including the\nexpanded description of the stressor task, the inclusion of HRV data as an additional physiological measure,\nand the discussion of key task-related limitations. We also appreciate the constructive identification of\nremaining areas that require clarification and additional detail and address these issues below.\nDespite the authors' points of clarification and justification, this reviewer is not entirely convinced regarding\nthe implementation and interpretation of results from the stress manipulation. A stressor of 2 minutes is very\nshort; while Dedovic et al (2005) had some subjects who had runs as short as 2 minutes, this was the\nminimum – most had runs >2 minutes. The bottom line is that this brief stressor has not been validated\nagainst the gold standard physiological metric: cortisol. This reviewer sees this point as a difference of\nopinion; while the authors may view HPA-axis measures as “complementary”, a review of the stress literature\nwould suggest that measures of cortisol are required to validate laboratory stressors (see, for example,\nDedovic et al, 2005, among others). As such, the limitation in the discussion should be rephrased from\n“Finally, the stress induction period was relatively brief (two minutes), and although we observed both\nsubjective and physiological signs of stress, neuroendocrine measures were not collected to provide\nadditional confirmation of hypothalamic-pituitary-adrenal-axis engagement” to “Finally, the stress induction\nperiod was relatively brief (two minutes), and although we observed both subjective and physiological signs\nof stress, this stress induction has not been validated using neuroendocrine measures (i.e., cortisol, ACTH),\nthus confirmation of hypothalamic-pituitary-adrenal-axis engagement is still needed.”\nResponse:\nWe agree with the Reviewer’s point and have revised the Discussion to reflect this more cautious\ninterpretation, incorporating the suggested wording regarding the lack of neuroendocrine validation of the\nstress induction.\nMS, p. 10:\n“Third, the stress induction period was relatively brief (two minutes), and although we observed both\nsubjective and physiological signs of stress, our stress induction task has not been validated against changes\nin cortisol or adrenocorticotropic hormone, which could provide confirmation of hypothalamic–pituitary–\nadrenal axis engagement70,71.\nThe title now uses the phrase “during high-stress” – “high-stress” relative to what? Given the results of the\nobjective ratings (mean of 5 on a 7-point scale) and short duration of the “stressor” (2 minutes), this\nterminology is misleading. I would recommend removing the term “high-stress” from the title and rephrasing\nto something like “Bed nucleus of the stria terminalis connectivity during food cue and taste processing in\nresponse to an acute laboratory stressor”. Additionally, the authors should acknowledge as a limitation that\nalthough subjective ratings between “low” and “high” stress differed, the ‘high” stress evoked only moderate\nstress levels on the 7-point rating scale.\nResponse:\nWe thank the Reviewer for this suggestion. We agree that “high-stress” may overstate the intensity of the\nmanipulation and have revised the title to use more neutral wording (“under stress”). Likewise, we have opted\nnot to use qualitative descriptors such as “moderate,” as these imply categorical thresholds not established\nin this study. For the sake of clarity, we have retained the terms “low-stress” and “high-stress” to describe the\nstudy’s conditions, without making claims about the intensity of stress experienced by participants.\nTitle:\n“Bed nucleus of the stria terminalis connectivity during food cue and taste processing under stress”\nRegarding the EDA data, text in the main document and supplementary figure are still misleading with regard\nto differences between “high” and “low” stress.\n- Page 5 of main document, paragraph starting “Physiological stress was assessed…”. Please\nrephrase the second sentence to “EDA values were not significantly different after high-stress…”\n(instead of “not significantly higher”).\n- Figure S1 caption: There is still commentary here regarding a stress-induced shift in EDA response\nwhich is misleading. Please revise lines 222-226 to \"The smoothed density plots (blue) show the\nestimated distribution of values in each condition. There was no difference between low- and high-\nstress conditions (t(27) = -0.31, p = 0.76, n = 28, mean difference = -0.006 [-0.05, 0.04]).\"\nResponse:\nWe thank the Reviewer for this suggestion. We have revised the text and figure caption as recommended to\nmore accurately reflect the absence of significant differences in electrodermal activity between the low- and\nhigh-stress conditions.\nIt is commendable that the authors utilized the pulse oximetry data to examine HRV as a physiological proxy\nfor stress. However, the manner of analysis and interpretation of the results remains unclear.\n- Why not examine RMSSD over each discrete task as a whole (i.e., the entire low-stress period; the\nentire taste task period following the low-stress induction, the entire high-stress induction; the entire\ntaste task period following the high-stress induction)? This should be done as an overall validity check,\nsimilar to the EDA data. That is, the authors should examine whether RMSSD differed across the\nwhole “high stress” run compared to across the whole “low stress” run.\nResponse:\nWe appreciate the reviewer’s positive feedback regarding the inclusion of pulse oximetry data to examine\nHRV. As the pulse trace could only be collected and synchronized during the acquisition of functional images\n(i.e., the task itself, not the stress induction phases), it is only possible to compare HRV measures from the\nbeverage task runs. We now explicitly state this in the Methods section:\nMS, p. 17:\n“Pulse oximetry was only collected and synchronized during functional image acquisition, meaning HRV\nestimates could only be derived from the beverage task runs.”\nOverall RMSSD values did not differ significantly between the low- and high-stress runs (t = −0.60, p =\n(43)\n0.55), indicating no detectable change in short-term heart-rate variability when averaged across the entire\nrun. However, examining overall RMSSD is known to mask meaningful stress responses as it collapses over\nperiods of varying physiological load and recovery within each task run. Stress reactivity typically unfolds\ndynamically, and heightened arousal may occur in specific segments of the task and later return toward\nbaseline. Analysing RMSSD in shorter, moving windows provides a more sensitive view of these temporal\npatterns, capturing transient changes can be lost when averaging across the full run. This windows approach\nis outlined in greater detail in our following response and this non-significant result for the whole run is now\nreported:\nMS, p. 17:\n“As stress responses unfold dynamically within a run, collapsing RMSSD across the entire task can obscure\nperiods of heightened arousal and subsequent recovery, whereas computing RMSSD in shorter sliding\nwindows provides the sensitivity needed to detect these transient physiological shifts. Indeed, mean RMSSD\ndid not differ when examining the entirety of low- and high-stress runs (t = −0.60, p = 0.55), whereas\n(43)\nsignificant differences did emerge using a sliding-windows approach.”\n- Instead, authors report that “Mean RMSSD HRV values during the beginning of the high-stress run\nwere lower than during the beginning of the low stress run (t(43) = -2.25, p = 0.029), indicating an\nelevated physiological stress response.” It is entirely unclear how this last point is concluded. It is also\nclaimed that the “high-stress condition was associated with sustained reductions in HRV” – sustained\nhow? Until when? It appears that some “beginning” of the run was defined for each version of the\nstress task, but how was “the beginning of the high-stress” run operationalized? And what does it\nmean if “the beginning” is different?\nResponse:\nWe thank the reviewer for raising these points. As detailed in the manuscript: The “beginning of the run”\nrefers to the first 20 seconds of each run, consistent with prior studies that assess acute autonomic responses\nto stress onset. RMSSD was computed within these 20-second windows to capture changes following task\ninitiation. The finding that RMSSD values were lower during the beginning of the high-stress run compared\nto the low-stress run (t(43)=−2.25, p=0.029) indicates a reduced short-term heart-rate variability at stress\nonset, reflecting an acute physiological stress response.\nTo address Reviewer’s 1 question in the first revision regarding sustained effects, we analysed the entire\ntime course of RMSSD using a sliding-window approach (10-second windows, 2-second steps) across the\nentire task. This dynamic analysis revealed a significant Condition × Time interaction (F(2,86)=3.20,\np=0.045), showing that RMSSD evolved differently across time between stress conditions. RMSSD values\nwere extracted continuously across the task period and subsequently segmented into early, middle, and late\ntask periods (first, middle, and final 20 seconds) to evaluate trends across time. The within-participants\nrepeated-measures ANOVA showed a significant Condition × Time interaction F(2, 86) = 3.20, p = 0.045,\nindicating that RMSSD evolved differently over time between the low- and high-stress conditions. Additional\npaired t-test comparisons showed that RMSSD increased significantly more from early to middle in the LS\ncondition compared to the high-stress condition, t(43) = 2.25, p=0.029. Similarly, from early to late, RMSSD\nincreased in the low-stress condition but decreased in the high-stress condition, t(43)=2.03, p=0.048.\nAs stated above, this approach examining RMSSD dynamics over successive windows provides a more\nsensitive and temporally precise measure of stress reactivity than averaging across an entire run. See:\nLi K, Rüdiger H, Ziemssen T. Spectral Analysis of Heart Rate Variability: Time Window Matters. Front Neurol.\n&\nPereira T, Almeida PR, Cunha JPS, Aguiar A. Heart rate variability metrics for fine-grained stress level\nassessment. Comput Methods Programs Biomed. 2017 Sep;148:71-80. doi: 10.1016/j.cmpb.2017.06.018.\nPMID: 28774440.\n- Next, the authors used a moving window analysis approach, citing a single paper that utilized this\napproach, but do not present a coherent reason for using this approach vs. RMSSD over the whole\n7-minute task. What is the validity of examining early, middle, and late task periods? It is unclear what\nit means that HRV was low in the early part of the task and higher later, given how short these “early”,\n“middle” and “late” phases are presumed to be.\nResponse:\nAveraging RMSSD over the full 7-minute task assumes stationarity and obscures transient autonomic shifts\nthat occur at onset and during recovery. The d=additional references supporting this approach have now\nbeen added to the manuscript. Short-time, moving-window analyses are recommended precisely because\nHRV (including RMSSD) is time-varying; window length materially affects what autonomic dynamics are\ncaptured, and shorter segments have been found to be better suited to track rapid changes than whole-\nrecording averages. Segmenting the task into early, middle, and late periods operationalises these dynamics:\nan “early” window indexes immediate reaction after stress onset, whereas later windows capture recovery or\nsustained effects. This position is now further supported in the methods section of the manuscript:\nMS, p. 18:\n“As stress responses unfold dynamically within a run, collapsing RMSSD across the entire task can\nobscure periods of heightened arousal and subsequent recovery, whereas computing RMSSD in shorter\nsliding windows provides the sensitivity needed to detect these transient physiological shifts. Indeed, mean\nRMSSD did not differ when examining the entirety of low- and high-stress runs (t = −0.60, p =0.55),\n(43)\nwhereas significant differences did emerge using a sliding-windows approach.\nRMSSD values were extracted continuously across the task period and subsequently segmented into early,\nmiddle, and late phases, defined as the first, middle, and final 20 seconds of the task, respectively. These\nvalues were entered into a within-participants repeated-measures ANOVA with Condition (low-stress, high-\nstress) and Time (early, middle, late) as factors. This design follows recommendations to model HRV as a\ntime-varying signal, allowing condition effects to be tested while capturing the distinct phases of autonomic\nreactivity within each run110,111.”\n- Relatedly, there are no figures to provide visualization of these complex data, which is strange given\nall of the other figures provided in support of other data in the manuscript. It is unclear how the data\nwere divided between the 2-minute stress blocks and the taste paradigm, and how the\nearly/middle/late task periods were defined. Showing plots of these data against the background of\nthe task would be helpful.\nResponse:\nWe thank the reviewer for this suggestion and agree that it would be helpful for readers to be provided with\na depiction of our HRV analyses. These are now included in Figure S2. This figure also includes further\ninformation on how data were collected during the fMRI taste task runs and how the early/middle/late task\nperiods were defined.\nSI, p. 18:\nTime course and phase differences in heart rate variability during low- and high-stress runs of the fMRI\nbeverage task. (a) Root mean square of successive differences (RMSSD) was extracted from pulse signals\nusing the PhysIO toolbox and calculated over 20-second windows stepped every 2 seconds. Shaded (grey)\nregions indicate early, middle, and late 20-second windows of the task. A repeated-measures ANOVA with\nfactors Condition (low-stress, high-stress) and Time (early, middle, late) revealed a significant Condition ×\nTime interaction, F = 3.20, p = 0.045, indicating differential RMSSD trajectories across stress conditions.\n(2, 86)\n(b) RMSSD increased significantly from early to middle and early to late in the low-stress run but declined in\nthe high-stress run (early → middle: t = 2.25, p = 0.029; early → late: t = 2.03, p = 0.048). Asterisks\n(43) (43)\ndenote significant between-condition differences.\nRegarding the SAGER guidelines, I might have missed it, but I do not see adherence to the requirement\n“Data should be reported disaggregated for sex and gender where this information has been collected.” I\nbelieve that data on sex was collected (Table S1), so disaggregated data should be presented. This is\nparticularly critical given the focus of the current manuscript on stress responsivity and connectivity in\nresponse to taste stimuli in several brain regions, which in prior literature have shown significant effects of\nsex (see references below). Notably, the BNST is highly sexually dimorphic and has been shown to have\nconnectivity patterns which differ by sex (see below). These extend further to include differences (in\nfemales) according to regular menstrual cycle variation (endogenous estradiol level) and use of oral\ncontraception and hormonal IUDs (see references below). Further, I do not see report of whether there was\nsufficient power to conduct specific sex-based analyses (i.e., not disaggregated, but directly compared),\nand if power is adequate, results of those analyses should be included, regardless of the outcome\n(significant differences or not). Finaly, it is unfortunate that gender data was not collected; this should be\nnoted as a limitation.\n- Bürger Z, Kordowich C, Kübbeler J, Müllerschön C, Kimmig AS, Su M, Lämmerhofer M, Sacher J,\nHenes M, Comasco E, Derntl B, Kogler L. Subjective, behavioural and physiological correlates of\nstress in women using hormonal contraceptives. Br J Psychiatry. 2025 Jun;226(6):392-400. doi:\n10.1192/bjp.2025.7. Epub 2025 Jun 13. PMID: 40511505; PMCID: PMC12257283.\n- Mengelkoch S, Gassen J, Slavich GM, Hill SE. Hormonal contraceptive use is associated with\ndifferences in women's inflammatory and psychological reactivity to an acute social stressor. Brain\nBehav Immun. 2024 Jan;115:747-757. doi: 10.1016/j.bbi.2023.10.033. Epub 2023 Oct 31. PMID:\n37914104; PMCID: PMC11216059.\n- Gervasio J, Zheng S, Skrotzki C, Pachete A. The effect of oral contraceptive use on cortisol reactivity\nto the Trier Social Stress Test: A meta-analysis. Psychoneuroendocrinology. 2022 Feb;136:105626.\ndoi: 10.1016/j.psyneuen.2021.105626. Epub 2021 Dec 11. PMID: 34922094.\n- Lovallo WR, Cohoon AJ, Acheson A, Vincent AS, Sorocco KH. Cortisol stress reactivity in women,\ndiurnal variations, and hormonal contraceptives: studies from the Family Health Patterns Project.\nStress. 2019 Jul;22(4):421-427. doi: 10.1080/10253890.2019.1581760. Epub 2019 Mar 21. PMID:\n30896268; PMCID: PMC6525075.\n- Aleknaviciute J, Tulen JHM, De Rijke YB, Bouwkamp CG, van der Kroeg M, Timmermans M, Wester\nVL, Bergink V, Hoogendijk WJG, Tiemeier H, van Rossum EFC, Kooiman CG, Kushner SA. The\nlevonorgestrel-releasing intrauterine device potentiates stress reactivity. Psychoneuroendocrinology.\n- Liu JJW, Ein N, Peck K, Huang V, Pruessner JC, Vickers K. Sex differences in salivary cortisol\nreactivity to the Trier Social Stress Test (TSST): A meta-analysis. Psychoneuroendocrinology. 2017\nAug;82:26-37. doi: 10.1016/j.psyneuen.2017.04.007. Epub 2017 Apr 24. PMID: 28486178.\n- Bürger Z, Müller VI, Hoffstaedter F, Habel U, Gur RC, Windischberger C, Moser E, Derntl B, Kogler\nL. Stressor-Specific Sex Differences in Amygdala-Frontal Cortex Networks. J Clin Med. 2023 Jan\n22;12(3):865. doi: 10.3390/jcm12030865. PMID: 36769521; PMCID: PMC9918214.\n- Brugnera A, Zarbo C, Tarvainen MP, Marchettini P, Adorni R, Compare A. Heart rate variability during\nacute psychosocial stress: A randomized cross-over trial of verbal and non-verbal laboratory\nstressors. Int J Psychophysiol. 2018 May;127:17-25. doi: 10.1016/j.ijpsycho.2018.02.016. Epub 2018\nMar 1. PMID: 29501671.\n- Cohen JE, Holsen LM, Ironside M, Moser AD, Duda JM, Null KE, Perlo S, Richards CE, Nascimento\nNF, Du F, Zuo C, Misra M, Pizzagalli DA, Goldstein JM. Neural response to stress differs by sex in\nyoung adulthood. Psychiatry Res Neuroimaging. 2023 Jul;332:111646. doi:\n10.1016/j.pscychresns.2023.111646. Epub 2023 Apr 25. PMID: 37146439; PMCID: PMC10247431.\n- Albert K, Pruessner J, Newhouse P. Estradiol levels modulate brain activity and negative responses\nto psychosocial stress across the menstrual cycle. Psychoneuroendocrinology. 2015 Sep;59:14-24.\ndoi: 10.1016/j.psyneuen.2015.04.022. Epub 2015 May 7. PMID: 26123902; PMCID: PMC4492530.\n- Ponticorvo S, Prinster A, Cantone E, Di Salle F, Esposito F, Canna A. Sex differences in the taste-\nevoked functional connectivity network. Chem Senses. 2022 Jan 1;47:bjac015. doi:\n10.1093/chemse/bjac015. PMID: 35749468.\n- Haase L, Green E, Murphy C. Males and females show differential brain activation to taste when\nhungry and sated in gustatory and reward areas. Appetite. 2011 Oct;57(2):421-34. doi:\n10.1016/j.appet.2011.06.009. Epub 2011 Jun 21. PMID: 21718731; PMCID: PMC3516392.\n- Smeets PA, de Graaf C, Stafleu A, van Osch MJ, Nievelstein RA, van der Grond J. Effect of satiety\non brain activation during chocolate tasting in men and women. Am J Clin Nutr. 2006 Jun;83(6):1297-\n305. doi: 10.1093/ajcn/83.6.1297. PMID: 16762940.\n- Flook EA, Feola B, Avery SN, Winder DG, Woodward ND, Heckers S, Blackford JU. BNST-insula\nstructural connectivity in humans. Neuroimage. 2020 Apr 15;210:116555. doi:\n10.1016/j.neuroimage.2020.116555. Epub 2020 Jan 16. PMID: 31954845; PMCID: PMC7089680.\n- Avery SN, Clauss JA, Winder DG, Woodward N, Heckers S, Blackford JU. BNST neurocircuitry in\nhumans. Neuroimage. 2014 May 1;91:311-23. doi: 10.1016/j.neuroimage.2014.01.017. Epub 2014\nJan 18. PMID: 24444996; PMCID: PMC4214684.\n- Allen LS, Gorski RA. Sex difference in the bed nucleus of the stria terminalis of the human brain. J.\nComp. Neurol. 1990;302:697–706. doi: 10.1002/cne.903020402.\n- Chung WCJ, De Vries GJ, Swaab DF. Sexual differentiation of the bed nucleus of the stria terminalis\nin humans may extend into adulthood. J. Neurosci. 2002;22:1027–1033. doi:\n10.1523/JNEUROSCI.22-03-01027.2002.\nResponse:\nWe thank the Reviewer for this point. We have added information on contraceptive use to Table S5 and note\nthat self-reported sex, but not gender, was collected.\nWe agree that prior work has reported sex-related differences in BNST size and connectivity, and we have\nrevised the Discussion to highlight this literature (see below). However, we note that whether such differences\nconstitute a “sexual dimorphism” is contentious (Eliot et al., 2021). In line with the Nature Communications\nreporting guidelines, we have not conducted post hoc sex-based analyses, as we do not believe the study\ndesign is sufficient to enable meaningful conclusions vis-à-vis sex-based differences. Specifically, we use a\nBayesian modelling framework for which conventional null-hypothesis-significance-testing power calculations\nare not applicable. Moreover, the relationship between sample size and reliable detection of posterior group\ndifferences (for example, thresholding at 0.75 posterior probability) has not been systematically established\n(cf. Ma et al., 2024).\nReferences:\nMa, L. et al. Effect of scanning duration and sample size on reliability in resting state fMRI dynamic causal\nmodeling analysis. NeuroImage 292, 120604 (2024).\nEliot, L., Ahmed, A., Khan, H. & Patel, J. Dump the “dimorphism”: Comprehensive synthesis of human brain\nstudies reveals few male-female differences beyond size. Neuroscience & Biobehavioral Reviews 125, 667–\n697 (2021).\nMS, p. 10:\n“Finally, although the BNST differs in size and connectivity between sexes and prior work highlights sex and\nhormonal influences on stress responsivity and gustatory processing 77–80, our study was not designed or\npowered to test these effects directly.”\nSI, p. 19:\nTable S5. Participant characteristics\nn = 48\nSex, female 50.0% (24)\nAge, years Mean = 26.2 ± 5.3\nBMI, kg/m2 Mean = 22.4 ± 3.8\nEthnicity\nAustralian 16.7% (8)\nBritish or European 18.8% (9)\nAsian 68.8% (33)\nMiddle Eastern 2.1% (1)\nAfrican 2.1% (1)\nNorth American 2.1% (1)\nCentral or Southern American 6.3% (3)\nEducation level\nYear 12 or equivalent 20.8% (10)\nBachelor’s degree 29.2% (14)\nPostgraduate degree 50.0% (24)\nHandedness, right 93.8% (45)\nDays since last period Mean = 16.4 ± 9.9\nContraceptives\nOral contraceptives 6.25% (3)\nIUD 0% (0)\nOther 2.08% (1)\nBMI, body mass index; IUD, intrauterine device. Data presented as % (n) or mean\n± standard deviation (SD). Days since last period are reported only for participants\nwho menstruate (40%, n = 19 of the total sample). Sex was self-reported, data on\ngender was not collected. \nResponses to the Reviewers\nTitle:\nBed nucleus of the stria terminalis connectivity during food cue and taste processing under stress\nManuscript: NCOMMS-25-24841B\nFor convenience, all comments are provided in boxed format with responses below. Quoted text from the\nmanuscript, supporting information and reporting summary appears in blue.\nReviewer 1\nI thank the authors for providing additional HRV measures as well as adjusting the language to better reflect\nthe post hoc nature of the analyses. Although I agree that the changes in HRV parameters support the\ninterpretation, changes in HR are, again, not reported.\nReply:\nWe thank the Reviewer for their thoughtful input throughout the review process. For completeness, we have\nnow included mean heart rate alongside the expanded HRV analyses. These analyses have been added to\nthe GitHub repository.\nMS, p. 6:\nNo significant differences were found in mean heart rate between the stress conditions: t = -0.37, p =\n(47)\n0.710.\nReviewer 2\nThe authors present an updated manuscript investigating BNST effective connectivity during taste under\nstress. The addition of the “over-complete” model was helpful to the reviewer and will be further helpful to the\nreaders. Some clarifications remain to be addressed:\n- Please briefly describe the results of the additional over-complete models for the readers (p. 7).\n- Provide figures for the over-complete models in the Supplemental file. This will help readers further compare\nbetween the primary and over-complete models.\n- Figure 2 left side panels are difficult to read—axis and labels are missing (presumably obscured by the\nblack background).\nResponse:\nWe thank the Reviewer for these helpful suggestions. We have addressed all points as requested: (i) we now\nbriefly summarise the results of the over-complete models in the main text; (ii) we have added figures for the\nover-complete models to the Supplementary Information to facilitate comparison with the primary models;\nand (iii) we confirm that all axis ticks and labels are legible in Fig. 2, and have rotated the ‘cue’ and ‘taste’\npanel identifiers for greater legibility. The relevant changes are reproduced below.\nMS, p. 7:\nResults from additional over-complete models—specified identically to the main analyses except that the\nmodulatory inputs were coded to represent all four conditions simultaneously—are also reported in the\nSupporting Information (Table S5 & Figure S4). These analyses suggest that modulation is primarily stress-\ndriven, with minimal evidence for stimulus-specific (water versus chocolate milk) effects.\nSI, p. 20:\nFigure S4. Stress-dependent modulation of BNST-centred effective connectivity during cue and taste\nconditions\nSchematic summarising group-level modulatory effects on directed effective connectivity (Hz) estimated with\nDCM when all four task conditions were modelled simultaneously. (a) Cue. (b), Taste. Only connections\nshowing supra-threshold modulation are displayed (posterior probability, Pp ≥ 0.75). Numbers on edges\nindicate the magnitude of modulation (±Hz), where the positive sign corresponds to low-stress (LS) and the\nnegative sign to high-stress (HS), reflecting the orthogonal Stress regressor coded as ±0.5; thus, the same\nconnection is upregulated under LS and sign-reversed under HS. In both epochs, stress modulated reciprocal\ncoupling between BNST and NAcc/dmINS and BNST–OFC coupling, with generally larger effects during cues\nthan during taste. The orthogonal Stimulus regressor (water vs chocolate milk; ±0.5) did not survive the Pp ≥\n0.75 criterion in either model, yielding identical condition-specific estimates within stress level (LS-Water =\nLS-Milk; HS-Water = HS-Milk). Region abbreviations: BNST, bed nucleus of the stria terminalis; NAcc,\nnucleus accumbens; OFC, orbitofrontal cortex; dmINS, dorsal mid-insular cortex.\nMS, p. 20:\nFigure 2. Behavioural ratings and group-level whole-brain general linear model activation during food\ncue and taste processing. (a) Mean self-reported hunger and stress ratings during low-stress and high-\nstress inductions, paired t-tests, *p < .05, ****p < .0001 (low stress vs high stress). Error bars show standard\nerror of the mean. (b) Mean self-reported chocolate milk and water pleasantness ratings after low- and high-\nstress inductions, 2-way ANOVA (stress x beverage), main effect of beverage, ****p < .0001, post hoc multiple\ncomparisons, paired t-tests, *p < .05, **p < .01 (low-stress vs high-stress). (c) Time course showing mean\nself-reported chocolate milk and water pleasantness ratings during each beverage cue exposure, 3-way\nANOVA (stress x beverage x cue), post hoc multiple comparisons unpaired t-tests at each cue exposure, ***p\n< .001, ****p < .0001 (chocolate milk vs water during high stress); ^^^p < .001; ^^^^p < .0001 (chocolate milk\nvs water during low stress); n = 48 for all. Error bars show standard error of the mean. (d) Group-level GLM\nresults for cue presentation (CUE), n = 44. (e) Group-level GLM results for taste receipt (TASTE), n = 47.\nGroup-level activation t-maps were overlaid on a standard MNI space template. The t-values are represented\nwith the colour bar. Significant activation was thresholded at a voxel-wise pFDR < .05, KE ≥ 10. Voxels of\ninterest are highlighted with yellow-dotted circles, including the OFC, dmINS, BNST and NAcc. An inclusive\nbinary mask was used during the NAcc time series extraction from the MNI-standardized automated\nanatomical labelling 3 (AAL3) toolbox, shown in red, and BNST time series extraction (see Methods) shown\nin yellow. BNST, bed nucleus of the stria terminalis; dmINS, dorsal mid-insula; NAcc, nucleus accumbens;\nOFC, orbitofrontal cortex; HS, high stress; LS, low stress.\nReviewer 3\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the\nNature Communications initiative to facilitate training in peer review and to provide appropriate recognition\nfor Early Career Researchers who co-review manuscripts.\nReply:\nAgain, we thank the Reviewer for their input in the review process.\nReviewer 4\nThe authors have thoughtfully responded to all of my comments, and the revisions, especially adding clarity\nand figures on the RMSSD data, have strengthened the manuscript.\nMy only remaining concern is regarding the SAGER guidelines. As stated in the Nature Research Ethics\nguidelines for Sex, Gender (identity/Presentation), and sexual orientation, parts iii and v:\niii. data should be reported disaggregated for sex and gender where this information has been collected and\nconsent has been obtained for reporting and sharing individual-level data; disaggregated numbers for\nindividual experiments must be provided in the source data as appropriate whereas overall numbers may be\nprovided in the Nature Portfolio Reporting Summary.\nv. If no sex- and gender-based analyses have been performed, authors should justify reasons for lack of\nanalysis in the Nature Portfolio Reporting Summary.\nThe authors responded that “We thank the Reviewer for this point. We have added information on\ncontraceptive use to Table S5 and note that self-reported sex, but not gender, was collected. We agree that\nprior work has reported sex-related differences in BNST size and connectivity, and we have revised the\nDiscussion to highlight this literature (see below). However, we note that whether such differences constitute\na “sexual dimorphism” is contentious (Eliot et al., 2021). In line with the Nature Communications reporting\nguidelines, we have not conducted post hoc sex-based analyses, as we do not believe the study design is\nsufficient to enable meaningful conclusions vis-à-vis sex-based differences. Specifically, we use a Bayesian\nmodelling framework for which conventional null-hypothesis-significance-testing power calculations are not\napplicable. Moreover, the relationship between sample size and reliable detection of posterior group\ndifferences (for example, thresholding at 0.75 posterior probability) has not been systematically established\n(cf. Ma et al., 2024).”\nHowever, I am not sure that the authors followed the suggestion as intended. That is, the response still does\nnot adhere to part iii above, encouraging disaggregation of data by sex. That is, reporting of results for males\nand for females, at least in supplementary materials. Citing a single study that casts doubt on sexual\ndimorphism of the brain does not negate the overwhelming number of studies that support sexual dimorphism\nrelated to stress and the BNST (listed in the prior round of reviews).\nThere is also no compliance with part v, wherein the authors should have justified the lack of post-hoc\nanalyses by sex in the Reporting Summary. Instead, the Reporting Summary merely states that equal\nnumbers of males and females were included. I do not see a justification for the lack of analysis by sex in the\nReporting Summary.\nResponse:\nWe thank the Reviewer for raising these points, which we are happy to clarify. Adherence to guideline (iv)\nnecessarily precludes the post hoc sex- or gender-disaggregated analyses implied by (iii) in this study, given\ninsufficient power:\niv. if sex- and gender-based analyses have been performed a priori, results should be reported\nregardless of positive or negative outcome. Authors should refrain from conducting post hoc sex- and\ngender-based analysis if the study design is insufficient (for example, low sample size) to enable\nmeaningful conclusions.\nWe have now updated the relevant section of the Reporting Summary (see below) so that it complies with\nguideline (v). Furthermore, we have now provided self-reported sex information in the GitHub repository (and\nnote this in the data and code availability statement (see below). We thank the Reviewer for their helpful input\nthroughout the review process.\nReporting Summary:\nThe aim of the study was to characterise BNST connectivity underpinning food-cue presentation and receipt\nfollowing acute stress exposure in the general population. Equal numbers of male and female participants\nwere recruited (n = 24 per sex; sex assigned at birth, self-reported) and are reported in the manuscript and\nSupporting Information; all participants provided informed consent for data sharing. Disaggregated data by\nsex are available in sex_ID.mat in the GitHub repository, allowing for separate sex-specific summary\nstatistics.\nThe study was not designed or powered for sex- or gender-based inferential analyses; accordingly, no post\nhoc sex/gender-disaggregated effective-connectivity analyses were conducted. Evidence suggests that in\nthe more statistically efficient resting-state effective connectivity setting—where explained variance exceeds\nthat of task-based analyses—reliable parameter recovery generally requires sample sizes >40 (Almgren et\nal., 2018; Zeidman et al., 2019; Ma et al., 2024).\nMS, p. 17–18:\nData and code availability. Deidentified effective connectivity, behavioural data, self-reported sex and\nscripts used to generate the main results of the study are publicly available at:\nhttps://github.com/evaghreins/Gustometer_Stress_BNST_DCM.\nResponses to the Reviewers\nTitle:\nBed nucleus of the stria terminalis connectivity during food cue and taste processing under stress\nManuscript: NCOMMS-25-24841C\nFor convenience, all comments are provided in boxed format with responses below. Quoted text from the\nmanuscript, supporting information and reporting summary appears in blue.\nReviewer 2\nThe authors have addressed the previous concerns.\nResponse:\nWe thank the Reviewer for their input in the review process.\nReviewer 3\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the\nNature Communications initiative to facilitate training in peer review and to provide appropriate recognition\nfor Early Career Researchers who co-review manuscripts.\nReply:\nAgain, we thank the Reviewer for their input in the review process.\nReviewer 4\nAdding the sex_ID.mat file to the GitHub is not a transparent, accessible way to report data disaggregated\nby sex. To actually report in a transparent, accessible way to the audience, one would include data\ndisaggregated by sex in the Supplementary Materials.\nResponse:\nSex per subject ID is now reported in the Supporting Information (see Table S6, reproduced below), and we\nrefer the reader to the study repository for further disaggregation (see ‘Demographic_Behaviour_data.xlsx’).\nThis change has been noted in the Reporting Summary. We thank the Reviewer for their input.\nSI, p. 11:\nN = 48\nSex, female 50.0% (24)\nAge, years Mean = 26.2 ± 5.3\nBMI, kg/m2 Mean = 22.4 ± 3.8\nEthnicity\nAustralian 16.7% (8)\nBritish or European 18.8% (9)\nAsian 68.8% (33)\nMiddle Eastern 2.1% (1)\nAfrican 2.1% (1)\nNorth American 2.1% (1)\nCentral or Southern American 6.3% (3)\nEducation level\nYear 12 or equivalent 20.8% (10)\nBachelor’s degree 29.2% (14)\nPostgraduate degree 50.0% (24)\nHandedness, right 93.8% (45)\nDays since last period Mean = 16.4 ± 9.9\nContraceptives\nOral contraceptives 6.25% (3)\nIUD 0% (0)\nOther 2.08% (1)\nBMI, body mass index; IUD, intrauterine device. Data presented as percentage or mean ±\nstandard deviation. Days since last period are reported only for participants who menstruate\n(40%, n = 19 of the total sample). Sex was self-reported, data on gender was not collected.\nParticipant IDs by self-reported sex; male (n = 24): 100, 101, 106, 111, 116, 117, 118, 119,\n121, 125, 127, 130, 132, 133, 134, 137, 138, 139, 142, 147, 150, 151, 152, 154; female (n =\n24): 102, 103, 104, 108, 109, 110, 113, 114, 115, 120, 122, 123, 126, 128, 131, 135, 140, 143,\n145, 146, 148, 149, 153, 156. See repository for further disaggregation."}]

Decision letter