Participants
Sixty-nine right-handed male participants with normal or corrected-normal vision were recruited from the University of Electronic Science and Technology of China. The sample size was determined using G*Power v.3.1 46, 47 indicating that 62 participants were required to detect a significant effect on the behavioral level (f\(=\)0.25, α\(=\)0.05, β\(=\)0.80, 2\(\times\)2\(\times\)2 mixed ANOVA interaction effect). The present study focused on male subjects to reduce variance related to sex differences in the effects of OT 48, 49 and adhered to validated exclusion criteria (details see supplementary methods). The final sample size was n\(=\)67 (mean\(\pm\)SD, age\(=\)21.03\(\pm\)1.95 years) due to two data from two participants being excluded (withdraw due to personal reasons, technical MRI issues). All participants provided written informed consent, protocols were pre-registered on Clinical Trials.gov (clinicaltrials.gov/ NCT05892939), approved by the ethics committee at the University of Electronic Science and Technology of China (UESTC-1061423041725893) and in line with the latest Declaration of Helsinki.
Using a double-blind randomized, placebo-controlled, between-subject pharmacological fMRI design, participants were randomly assigned to receive either a single intranasal dose of oxytocin (OT, 24IU-Sichuan Defeng Pharmaceutical Co. Ltd, China) or a placebo (PLC, same ingredients i.e. sodium chloride and glycerin but without the peptide). The spray bottles were identical and dispensed by an independent researcher based on a computer-generated randomization sequence to ensure double-blinding. The investigators involved in data acquisition and analyses were blinded for group allocation. Following recent recommendations 50–52, the fMRI assessment began 45 minutes after treatment administration.
To control for pre-treatment between-group differences, participants completed a series of mood and mental health questionnaires (as detailed in Table 1). Participants were asked to guess which treatment they had received after the experiment (with a non-significant treatment guess χ2\(=\)0.39, p\(=\)0.53, confirming successful double-blinding). Additionally, participants were required to complete a surprise recognition task after the experiment that involved presenting stills of the video clips shown during fMRI and new clips to test the attentive processing of the stimuli. All participants achieved a re-recognition accuracy of over 75% and no significant difference in performance was observed between the two treatment groups (t(1,65) \(=\)-0.693, p\(=\)0.491).
Naturalistic fear induction paradigm.
The task began 45 minutes after the treatment intervention (see Fig. 1a). During the naturalistic fear induction paradigm, participants were instructed to watch a series of video clips attentively and rate their level of fear experience after each clip. The clips lasted 25 seconds and were followed by a 1–3 second jittered fixation cross to separate the stimuli from the rating period. During the 5-second period following each stimulus, participants were asked to report their level of subjective fear on a 9-point Likert scale, with 1 indicating no fear and 9 indicating very strong fear. This was followed by a jittered inter-trial interval of 9–11 seconds, during which a fixation cross was presented (see Fig. 1b). A total of 32 short video clips were presented over 2 runs, with each run consisting of 16 randomly presented clips, half of which were neutral and half of which were fear-inducing clips.
Given ongoing discussions about the social-specific effects of OT 53, 54, we carefully pre-selected a balanced set of stimuli that included social as well as non-social fearful or neutral stimuli, respectively. This led to a set of four stimulus categories (fear social clips-FS, fear non-social clips-FNS, neutral social clips-NS, neutral non-social clips-NNS, see Fig. 1c). For detailed characterization of the video clips, see Supplementary Methods. The video clips used for the naturalistic fear induction paradigm were selected from a larger database based on stimuli evaluation ratings in an independent sample (n = 20; 10 females; age = 22.65 ± 2.01). Key exclusion criteria for the video clips were reported in supplementary methods. Based on the results of the behavioral ratings (see Fig. S2), 8 complimentary video clips were selected for each of the four categories. The task was programmed in Python 3.7 using the PsychoPy package (version 2022.2.4).
Effects of OT on subjective fear experience.
To determine the effects of OT on subjective fear in naturalistic contexts and whether they vary as a function of social context, we computed a repeated measures ANOVA with the between-subject factor treatment (OT vs PLC) and the within-subject factors emotion (fear vs neutral) and social (social vs non-social) and the dependent variable fear ratings during fMRI. To initially ensure that the selected video clips were effective in inducing the intended emotions, repeated measures ANOVA (sex was used as the independent variable, while emotion and social were used as repeated variables) were conducted on five rating metrics (see supplementary methods). Analyses were conducted using GraphPad Prism 9.0 (GraphPad Software, Inc., La Jolla, CA), and corresponding effect sizes were computed using JASP 0.18.3.0 (JASP Team, 2019; jasp-stats.org). All ANOVAs and t-tests used false discovery rate (FDR) 55 recommended by Prism for multiple comparison correction.
MRI data acquisition and preprocessing.
MRI data were collected on a 3.0-T GE Discovery MR750 system (General Electric Medical System, Milwaukee, WI, USA) and preprocessed using FMRIPREP 21.0.056 (RRID: SCR_016216), a Nipype 1.6.1 based tool that integrates preprocessing routines from different software packages and SPM 12 (Statistical Parametric Mapping; http://www.fil.ion.ucl.ac.uk/spm/; Wellcome Trust Centre for Neuroimaging, detailed in supplementary methods).
Next, a voxel-wise general linear model (GLM) was conducted for each participant. Specifically, a first-level model was designed that included separate regressors for the four conditions and the rating period. In line with a previous study 28, we removed variance associated with the mean, linear and quadratic trends, the average signal within anatomically-derived CSF mask, the effects of motion estimated during the head-motion correction using an expanded set of 24 motion parameters (six realignment parameters, their squares, their derivatives, and their squared derivatives) and motion spikes (FMRIPREP default: FD\(>\)0.5mm or standardized DVARS\(>\)1.5)57.
Effects of OT on fear-related brain activity.
To examine the effects of OT on fear-related neural activity in naturalistic contexts, a GLM-based model was established in SPM 12. The model incorporated five regressors: FS, FNS, NS, NNS, and the rating period. We included the 24 head motion parameters as covariates of no interest. Key first-level contrasts of interest were contrasts that modeled fear in social and non-social contexts separately (i.e., [FS\(-\)NS] or [FNS\(-\)NNS]) as well as their interaction (i.e., [FS\(-\)NS] \(-\)[FNS\(-\)NNS]). In line with the behavioral analyses, the effects of OT were determined using an ANOVA with the factors emotion and social context, implemented in a partitioned error ANOVA approach that subjected the first level interaction contrast to a voxel-wise two-sample t-test with the factor treatment group on the second level. Based on our hypothesis and previous studies reporting that fear is strongly associated with amygdala, cingulate cortex (particularly middle cingulate cortex, MCC), insula, and ventromedial prefrontal cortex (vmPFC) and that oxytocin mediates its anxiolytic effects via these regions 14, 22, 27, 58–60, our analyses focused on atlas-derived bilateral masks for the amygdala, MCC, insula and vmPFC using separate small volume correction (SVC) and family-wise error (FWE) correction 61 applied at peak level.
To further disentangle complex interaction effects involving treatment, parameter estimates were extracted using 8mm radius spheres centered at the peak interaction coordinates using Marsbar (http://marsbar.sourceforge.net/). The extracted Blood Oxygenation Level Dependent (BOLD) signals were then subjected to appropriate post-hoc tests comparing group differences in each of the separate conditions.
Additionally, to determine associations between the behavioral (subjective fear) and neural effects of OT (both activation and functional connectivity between specific brain regions) we conducted correlation analysis using JASP.
Effects of OT on functional communication of the regions identified.
We further examined whether the regions showing altered regional activation following OT also changed their network-level communication via a general Psychophysiological interaction (gPPI) analysis. gPPI examines functional connectivity between brain regions that is contingent on a psychological context. We used the Generalized PPI toolbox 62 as implemented in SPM 12, to conduct gPPI. This method includes additional task regressors to reduce the likelihood that the functional connectivity estimates are driven simply by co-activation and has been widely used in similar studies 63, 64. Based on our results, OT enhanced activation in the left middle cingulate cortex (lMCC) under the three-way interaction (treatment\(\times\)emotion\(\times\)social), we therefore computed the functional communication of this region (using an 8mm radius which centered in peak Montreal Neurological Institute (MNI), x, y, z = -7, -3, 30) with other voxels in the whole brain under the two-way interaction of emotion and social, and a two-sample t-test was conducted to confirm the treatment effect of OT. Based on previous studies reporting effects of OT on the functional communication of the amygdala 48, 49, separate atlas-derived bilateral masks for the amygdala, MCC, insula and vmPFC were used for SVC analyses with FWE correction 61 applied at peak level.
Marsbar was used to extract information communications between other regions and the seed region to further disentangle complex interaction effects involving treatment. The extracted functional connectivity signals were then subjected to compare the treatment effect under four conditions by using two-sample t-tests separately.
Effects of OT on large-scale brain network functional connectivity.
Given that several studies reported the effects of OT on large-scale network organization of the brain 65, 66, we examined OT effects on the whole-brain network level. To this end, time-series of each region of interest (ROI) was extracted according to the template including 463 ROIs across the whole brain (see supplementary methods), resulting in BOLD time series of the size of r\(\times\)t (r was the number of ROIs and t was the BOLD time series length) for each subject per run. A matrix of size r\(\times\)r was obtained by calculating Pearson's correlation between the average BOLD time series of each pair of ROIs. Additionally, Fisher's z-transform was applied to improve the normality of the correlation coefficients. The brain regions can be divided into nine networks based on their functional divisions, and the neural activities of each network's regions were averaged to calculate network functional connectivity.
During each run, four types of video clips were presented in a pseudo-randomized order. The onset time of each video clip was shifted by 3 TRs to account for hemodynamic lag. To extract neural activity time series for a specific condition, we used the time points and durations specified for the appearance of different trials of the same type. This resulted in four r\(\times\)d matrices (where r is the number of ROIs and d is the BOLD time series length of each condition). Additionally, we obtained the r\(\times\)r matrix between ROIs and the n\(\times\)n matrix (where n is the number of networks) between brain networks for different conditions.
After performing Fisher-Z transformation on the network functional connectivity matrices for different conditions, matrices of Z-values corresponding to specific conditions (ZFS\ZFNS\ZNS\ZFNS) are obtained. It is important to note that r values are not additive or subtractive, but Z scores are summative. Matrix of Z-values corresponding to particular contrasts (e.g. ZFS-NS) can be computed by subtract the matrices of Z-values corresponding to different conditions (e.g. ZFS-ZNS)67, 68. Therefore, each participant in both groups would have their own subtracted Z-values matrix under specific contrasts (ZFS-NS\ZFNS-NNS). We then conducted two-sample t-tests on each pair of network-level functional connectivity to determine the treatment effect.
Effects of OT on whole-brain functional connectivity.
To further investigate the impact of OT on a finer spatial scale, we analyzed whole-brain functional connectivity using paired t-tests for key contrasts (FS\(-\)NNS and FNS\(-\)NNS) separately for the two groups. We applied FDR correction (pFDR−corrected\(<\)0.05) to identify the most relevant functional connections for each contrast. A Mantel test 69, 70 was performed on the FDR-corrected functional connectivity matrices of the two groups, under specific contrasts, to assess their similarity. Furthermore, to investigate treatment-specific effects on functional connectivity strength at the whole-brain level, we transformed the upper triangular data of the t-value matrices of the two groups under specific contrast into two (r\(\times\)(r\(-\)1)) / 2 by 1 feature vectors, and performed Welch’s t-test between them.