Regulatory approvals and registrations
Written informed consent was obtained from all participants in accordance with the Declaration of Helsinki and procedures established by the Washington University in Saint Louis Institutional Review Board. All participants were compensated for their time. All aspects of this study were approved by the Washington University School of Medicine (WUSOM) Internal Review Board, the Washington University Human Research Protection Office (WU HRPO), the Federal Drug Administration (IND no. 202002165) and the Missouri Drug Enforcement Agency (DEA) under a federal DEA schedule 1 research licence and registered with ClinicalTrials.gov identifier NCT04501653. Psilocybin was supplied by Usona Institute through Almac Clinical Services.
Study design
Healthy young adults (n = 7, 18–45 years) were enrolled between April 2021 and March 2023 in a randomized crossover precision functional brain mapping study at Washington University in Saint Louis (see Supplementary Methods for inclusion and exclusion criteria). The purpose of the study was to evaluate differences in individuallevel connectomics before, during and after psilocybin exposure. Participants underwent imaging during drug sessions (with MRI starting 1 h after drug ingestion) with 25 mg psilocybin or 40 mg MTP, as well as nondrug imaging sessions. Drug condition categories were (1) baseline, (2) drug 1 (MTP or psilocybin), (3) between, (4) drug 2 and (5) after. Randomization allocation was conducted using REDCap and generated by team members who prepared study materials including drug or placebo but otherwise had no contact with participants. A minimum of three nondrug imaging sessions were completed during each nondrug window: baseline, between and after drug sessions. The number of nondrug MRI sessions was dependent on availability of the participant, scanner and scanner support staff. Dosing day imaging sessions were conducted 60–180 min following drug administration during peak blood concentration^{98}. One participant (P2) was not able tolerate fMRI while on psilocybin, and had trouble staying awake on numerous fMRI visits after psilocybin and was thus excluded from analysis (except for data quality metrics in Extended Data Fig. 1).
MTP was selected as the active control condition to simulate the cardiovascular effects and physiological arousal (that is, controlling for dopaminergic effects) associated with psilocybin^{99}. Usona Institute, a US nonprofit medical research organization, provided good manufacturing practices for psilocybin.
Drug sessions were facilitated by two clinical research staff who completed an approved inperson or online facilitator training programme provided by Usona Institute, as part of the phase 2 study (ClinicalTrials.gov identifier NCT03866174). The role of the study facilitators was to build a therapeutic alliance with the participant throughout the study, prepare them for their drug dosing days and to observe and maintain participant safety during dosing day visits^{64}. The pair consisted of an experienced clinician (lead clinical facilitator) and a trainee (cofacilitator).
The predefined primary outcome measure was precision functional mapping (numerous visits, very long scans to produce individual connectomes) examining the effects of psilocybin on cortical and cortico subcortical brain networks that could explain its rapid and sustained behavioural effects. Predefined secondary outcome measures included (1) assessment of hemodynamic response to evaluate how 5HT_{2A} receptor agonism by psychedelics may alter neurovascular coupling, (2) assessment of acute psychological effects of psilocybin using the MEQ30 score (Supplementary Methods) and (3) assessment of personality change using the International Personality Item PoolFiveFactor Model^{100}. Changes in pulse rate and respiratory rate during psilocybin and placebo were later added as secondary outcome measures and personality change was abandoned because it was clear that we would not be powered to detect personality change.
Replication protocol
Participants were invited to return 6–12 months after completing the initial crossover study for a replication protocol. This included 1–2 baseline fMRIs, a psilocybin session (identical to the initial session, except for lack of blinding) and 1–2 ‘after’ sessions within 4 days of the dose.
Participants
Healthy adults aged 18–45 years were recruited by campuswide advertisement and colleague referral. Participants (n = 7) were enrolled from March 2021 to May 2023. Participants were required to have had at least one previous lifetime psychedelic exposure (for example, psilocybin, mescaline, ayahuasca, LSD), but no psychedelics exposure within the past 6 months. Individuals with psychiatric illness (depression, psychosis or addiction) based on the DSM5 were excluded. Demographics and data summary details are provided in Supplementary Table 1. One of the authors (N.U.F.D.) was a study participant.
MRI
Participants were scanned roughly every other day over the course of the experiment (Extended Data Fig. 1). Imaging was performed at a consistent time of day to minimize diurnal effects in FC^{101}. Neuroimaging was performed on a Siemens Prisma scanner (Siemens) in the neuroimaging laboratories at the Washington University Medical Center.
Structural scans (T1w and T2w) were acquired for each participant at 0.9 mm isotropic resolution, with realtime motion correction. Structural scans from different sessions were averaged together for the purposes of Freesurfer segmentation and nonlinear atlas registrations.
To capture highresolution images of blood oxygenation leveldependent (BOLD) signal, we used an echoplanar imaging sequence^{102} with 2 mm isotropic voxels, multiband 6, multiecho 5 (times to echo: 14.20, 38.93, 63.66, 88.39, 113.12 ms)^{103}, repetition or relaxation time: 1,761 ms, flip angle of 68° and inplane acceleration^{104} (IPAT or grappa) of 2. This sequence acquired 72 axial slices (144 mm coverage). Each resting scan included 510 frames (lasting 15 min, 49 s) as well as three frames at the end used to provide estimate electronic noise.
Every session included two 15min restingstate fMRI (rsfMRI) scans, during which participants were instructed to hold still and look at a white fixation crosshair presented on a black background. Head motion was tracked in real time using Framewise Integrated Realtime MRI Monitoring software (FIRMM)^{105}. An eyetracking camera (EyeLink) was used to monitor participants for drowsiness.
Perceptual (matching) fMRI task
Participants also completed a previously validated eventrelated fMRI task. This was a suprathreshold auditory–visual matching task in which participants were presented with a naturalistic visual image (duration 500 ms) and coincident spoken English phrase, and were asked to respond with a button press to indicate whether the image and phrase were ‘congruent’ (for example, an image of a beach and the spoken word ‘beach’) or ‘incongruent’. Both accuracy and response time of button presses were recorded. Each trial was followed by a jittered interstimulus interval optimized for eventrelated designs. In a subset of imaging sessions, two task fMRI scans were completed following the two resting scans. Task fMRI scans used the same sequence used in resting fMRI, included 48 trials (24 congruent, 24 incongruent) and lasted a total of 410 s. In analyses, high motion frames were censored^{106} and the two task scans were concatenated to better match the length of the rsfMRI scans. Note the stimulus order in the two trials did not vary across session. The order of rest and task scans was not counterbalanced across sessions to avoid concern that task scans may influence subsequent rest scans.
Resting fMRI processing and restingstate network definition
Resting fMRI data were preprocessed using an inhouse processing pipeline. In brief, this included removal of thermal noise using NORDIC denoising^{107,108,109}, correction for slice timing and field distortions, alignment, optimal combination of many echoes by weighted summation^{110}, normalization, nonlinear registration, bandpass filtering and scrubbing at a movement threshold of 0.3 mm to remove reduce the influence of confounds^{111}. Tissuebased regressors were computed in volume (white matter, ventricles, extraaxial cerebrospinal fluid)^{112} and applied following projection to surface. Taskbased regressors were only applied when indicated. Details on rsfMRI preprocessing are provided in Supplementary Methods. Visualizations of motion, physiological traces and signal across the brain (‘grayplots’) before and after processing^{113} are provided in Supplementary Video 1.
Surface generation and brain areal parcellation
Surface generation and processing of functional data followed similar procedures to Glasser et al.^{114}. To compare FC and restingstate networks across participants, we used a groupbased surface parcellation and community assignments generated previously^{62}.
For subcortical regions, we used a set of regions of interest^{115} generated to achieve full coverage and optimal region homogeneity. A subcortical limbic network was defined on the basis of neuroanatomy: amygdala, anteromedial thalamus, nucleus accumbens, anterior hippocampus and posterior hippocampus^{116,117}. These regions were expanded to cover anatomical structures (for example, anterior hippocampus)^{31}.
To generate regionwise connectivity matrices, time courses of all surface vertices or subcortical voxels within a region were averaged. FC was then computed between each region timeseries using a bivariate correlation and then Fisher ztransformed for group comparison.
Individualized network and brain area mapping
We identified canonical largescale networks using the individualspecific network matching approach described previously^{43,44,62}. In brief, cortical surface and subcortical volume assignments were derived using the graphtheorybased Infomap algorithm^{118}. In this approach, we calculated the correlation matrix from all cortical vertices and subcortical voxels, concatenated across all a participant’s scans. Correlations between vertices within 30 mm of each other were set to zero. The Infomap algorithm was applied to each participant’s correlation matrix thresholded at a range of edge densities spanning from 0.01 to 2%. At each threshold, the algorithm returned community identities for each vertex and voxel. Communities were labelled by matching them at each threshold to a set of independent group average networks described previously^{62}. In each individual and in the average, a ‘consensus’ network assignment was derived by collapsing assignments across thresholds, giving each node the assignment it had at the sparsest possible threshold at which it was successfully assigned to one of the known group networks. See Extended Data Fig. 4 and Supplementary Fig. 1 for individual and group mode assignments, respectively. The following networks were included: the association networks including the DMN, frontoparietal, dorsal attention, parietal memory, ventral attention, actionmode, salience and context networks; and the primary networks including the visual, somatomotor, somatomotor face and auditory networks.
To compute local (areal) desynchronization, we also defined brain areas at the individual level using a previously described areal parcellation approach^{39}. In brief, for each participant, vertexwise FC was averaged across all sessions to generate a dense connectome. Then, abrupt transitions in FC values across neighbouring vertices were used to identify boundaries between distinct functional areas.
LME model
To take advantage of the multilevel precision functional mapping study design, a LME model was used. Every scan was labelled on the following dimensions: participant identity (ID), MRI visit, task (task or rest), drug condition (prepsilocybin, psilocybin, MTP, postpsilocybin) and head motion (average framewise displacement). The rsfMRI metrics (described below) were set as the dependent variable, drug (drug condition), task, framewise displacement (motion) and drug × task were defined as fixed effects, and participant ID and MRI session were random effects.
Let y_{ij} be the rsfMRI metric (for example, FC change score at a given vertex) for the jth observation (15 min fMRI scan) within the ith participant. The LME model can be written as:
$$\begin{array}{l}{y}_{ij}={\beta }_{0}+{\beta }_{{\rm{d}}{\rm{r}}{\rm{u}}{\rm{g}}}\cdot {{\rm{d}}{\rm{r}}{\rm{u}}{\rm{g}}}_{ij}+{\beta }_{{\rm{F}}{\rm{D}}}\cdot {{\rm{F}}{\rm{D}}}_{ij}+{\beta }_{{\rm{t}}{\rm{a}}{\rm{s}}{\rm{k}}}\cdot {{\rm{t}}{\rm{a}}{\rm{s}}{\rm{k}}}_{ij}\\ \,\,+{\beta }_{{\rm{task}} \mbox{} {\rm{by}} \mbox{} {\rm{drug}}}\cdot {{\rm{t}}{\rm{a}}{\rm{s}}{\rm{k}}}_{ij}\cdot {{\rm{d}}{\rm{r}}{\rm{u}}{\rm{g}}}_{ij}+{u}_{0i}+{v}_{0j}+{{\varepsilon }}_{ij}\end{array}$$
(1)

β_{0} is the intercept term.

β_{drug}, β_{FD}, β_{task} and β_{taskbydrug} are the coefficients for the fixed effects predictors.

drug_{ij}, frame displacement_{ij} (FD_{ij}) and task_{ij} are the values of the fixed effects predictors for the jth observation within the ith group.

u_{0i} represents the random intercept for the ith participant, accounting for individualspecific variability.

v_{0j} represents the random intercept for the jth observation within the ith participant, capturing scanspecific variability.

ε_{ij} is the error term representing unobserved random variation.
In MATLAB (Wilkinsonian notation), this model is expressed for every vertex Y(vertex) = fitlme(groupd, FC_Change(vertex) ~ drug + framewise displacement + task + taskbydrug + (1 SubID) + (1 session)).
To compensate for the implementations of this LME model on many rsfMRIrelated dependent variables, differences were highlighted when P < 0.001. All P values reported are not corrected for multiple comparisons.
Vertexwise FC change
FC change (‘distance’) was calculated at the vertex level to generate FC change maps and a LME model (equation (1)) was used in combination with wild bootstrapping^{119,120} and thresholdfree cluster enhancement (TFCE)^{95,121} to estimate P values for tstatistic maps resulting from the model (Figs. 1a–d and 4). Wild bootstrapping is an approach to permutation testing that was designed for models that are not independent and identically distributed, and are heteroscedastic.
First, a FC change map was generated for every scan by computing, for each vertex, the average distance between its FC seedmap and the FC seedmap for each of that participant’s baseline scans. As each participant had several baseline visits, FC change was computed for baseline scans by computing distance from all other baseline scans (excluding scans within the same visit). This provided a measure of daytoday variability. Second, the distance value was used as the dependent variable y_{ij} in the LME model to generate a tstatistic. Third, a wild bootstrapping procedure was implemented as follows. Several bootstrap samples (B = 1,000) were generated using the Rademacher procedure^{120}, in which the residuals were randomly inverted. Specifically, a Rademacher vector was generated by randomly assigning −1 or 1 values with equal probability to the residual of each observation. By elementwise multiplication of the original residuals with the Rademacher vector, bootstrap samples were created to capture the variability in the data.
For the observed tstatisticmap and each bootstrap sample, the TFCE algorithm was applied to enhance the sensitivity to clusters of significant voxels or regions while controlling for multiple comparisons. The value of the enhanced cluster statistic derived from the bootstrap samples was used to create a null distribution under the null hypothesis. By comparing the original observed cluster statistic with the null distribution, P values were derived to quantify the statistical significance of the observed effect. The P values were obtained on the basis of the proportion of bootstrap samples that produced a maximum cluster statistic exceeding the observed cluster statistic.
The combined approach of wild bootstrapping with the Rademacher procedure and TFCE provided the method to estimate P values for our multilevel (drug condition, participant, session, task) design. This methodology accounted for the complex correlation structure, effectively controlled for multiple comparisons and accommodated potential autocorrelation in the residuals through the Rademacher procedure. By incorporating these techniques, association with psilocybin and other conditions was reliably identified amid noise and spatial dependencies.
Wholebrain FC change
For analyses in Figs. 1e,g, 2 and 4a (bottom), Extended Data Fig. 3 and Supplementary Figs. 3, 4 and 6, distance calculations were computed on the FC matrix using ztransformed bivariate correlation of time courses from parcellated brain areas^{62}. The effects of daytoday, drug condition, task and framewise displacement and drug × task were directly examined by calculating the distance between functional network matrices generated from each scan. Rootmeansquared Euclidean distance was computed between the linearized upper triangles of the parcellated FC matrix between each pair of 15 min fMRI scans, creating a secondorder distance matrix (Extended Data Fig. 3). Subsequently, the average distance (reported as ‘wholebrain FC change’) was examined for FC matrices that were from the same individual within a single session, from the same individual across days (‘daytoday’), from the same participant between drug and baseline (for example, psilocybin), from the same individual but different tasks (‘task:rest’), from the same individual between highest motion scans and baseline (‘hi:lo motion’), from different individuals (‘between person’). In the ‘high head motion’ comparison (‘hi:lo motion’ in Supplementary Fig. 3), the two nondrug scans with the highest average framewise displacement were labelled and compared against all other baseline scans.
A LME model (equation (1)) and post hoc ttests were used to assess statistical differences between drug conditions. A related approach using ztransformed bivariate correlation (‘similarity’ rather than distance) was also taken and results were unchanged (Supplementary Fig. 3c).
Likelihood ratio test of participantspecific response
To test whether variability in participantspecific response to psilocybin was larger than would be expected by chance, we used a likelihood ratio test for variance of random slopes for a participantspecific response to psilocybin^{48}. The difference in log likelihood ratios was compared to a null distribution of 1 million draws from a mixture of chisquared distributions with degrees of freedom 1 and 2. We note that the likelihood ratio test of variance components is a nonstandard problem^{47} as the covariance matrix of the random effects is positive definite and the variances of random effects are nonnegative. Finally, the test statistic for the likelihood ratio in this LME model was compared against a 50/50 mixture of two independent chisquared distributions, each with one and two degrees of freedom, respectively.
Assessing subjective experience
Subjective experience was assessed for drug sessions using the MEQ30^{46} (Supplementary Methods). The MEQ30 is designed to capture the core domains of the subjective effects of psychedelics (as compared to the altered states of consciousness rating scales that more broadly assess effects of psychoactive drugs^{122}) and is related to the therapeutic benefits of psychedelics. We applied a LME model across all drug sessions, similar to the one described above, but with MEQ30 total score as the dependent variable. Wholebrain FC change and framewise displacement were modelled as fixed effects, and participant was modelled as a random effect. The same model was solved using FC change from every vertex to generate a vertexwise map of the FC change versus MEQ30.
Normalized FC change
The conditions above were compared by calculating normalized FC change scores using the following procedure: we (1) determined FC change for each condition compared to baseline as described above, (2) subtracted withinsession distance for all conditions (such that withinsession FC change was 0), (3) divided all conditions by daytoday distance (such that daytoday FC change was equal to 1). Thus, normalized wholebrain FC change values (for example, psilocybin versus base was 3.52) could be thought of as proportional to daytoday variability.
Datadriven MDS
We used a classical MDS approach to cluster parcellated connectomes across fMRI scans, as previously described^{38}. This datadriven approach was used to identify how different parameters (for example, task, drug, individual) affect similarity and/or distance between networks. MDS places data in multidimensional space on the basis of the dissimilarity (Euclidean distance) among data points, which in this case means a data point represents the linearized upper triangle of a FC matrix. Every matrix was entered into the classical MDS algorithm (implemented using MATLAB 2019, cmdscale.m). Many dimensions of the data were explored. The eigenvectors were multiplied by the original FC matrices to generate a matrix of eigenweights that corresponded to each dimension. These eigenweights were also applied to other rsfMRI psychedelics datasets to generate dimensions scores (section ‘Other datasets’).
Rotationbased null model (spin test) for network specificity
To assess network specificity of FC change values, we calculated average FC change of matched null networks consisting of randomly rotated networks with preserved size, shape and relative position to each other^{62,97}. To create matched random networks, we rotated each hemisphere of the original networks a random amount around the x, y and z axes on the spherical expansion of the cortical surface^{62}. This procedure randomly relocated each network while maintaining networks’ sizes, shapes and relative positions to each other. Random rotation followed by computation of networkaverage FC change score was repeated 1,000 times to generate null distributions of FC change scores. Vertices rotated into the medial wall were not included in the calculation. Actual psilocybin FC change was then compared to null rotation permutations to generate a P value for the 12 networks that were consistently present across every participant’s Infomap parcellation. For bar graph visualization (Fig. 1 and Supplementary Fig. 1b), networks with greater change (P < 0.05 based on null rotation permutations) are shown in their respective colour and other networks are shown in grey.
NGSC
We used an approach previously validated to assess spatial complexity (termed entropy) or neural signals^{61}. Temporal principal component analysis was conducted on the full BOLD dense timeseries, which yielded m principal components (m roughly 80 K surface vertices and subcortical voxels) and associated eigenvalues. The normalized eigenvalue of the ith principal component was calculated as
$${\lambda }_{i}^{{\prime} }=\frac{{\lambda }_{i}}{{\sum }_{i=1}^{m}{\lambda }_{i}^{{\prime} }}$$
(2)
where m is the number of principal components, and λ_{i} and λ′_{i} represent the eigenvalue and the normalized eigenvalue of the ith principal component, respectively. Last, the NGSC, defined as the normalized entropy of normalized eigenvalues, was computed using the equation:
$${\rm{NGSC}}=\frac{{\sum }_{i=1}^{m}{\lambda }_{i}^{{\prime} }\log {\lambda }_{i}^{i}\,}{\log m}$$
(3)
The NGSC computed above attains values from the interval 0 to 1. The lowest value NGSC = 0 would mean the brainwide BOLD signal consisted of exactly one principal component or spatial mode, and there is maximum global FC between all vertices. The highest value NGSC = 1 would mean the total data variance is uniformly distributed across all m principal components, and a maximum spatial complexity or a lowest FC is found.
NGSC was additionally calculated at the ‘parcel level’. To respect areal boundaries, this was done by first generating a set of individualspecific parcels in every participant (on all available resting fMRI sessions concatenated) using procedures described oreviously^{39,62}.
NGSC maps were compared to PETbased 5HT_{2A} receptor binding maps published in ref. ^{33}. Similarity was assessed by computing the bivariate correlation between NGSC values and 5HT_{2A} binding across 324 cortical parcels from the Gordon–Laumann parcellation.
Persistent effects analysis
To assess the persistent effects of psilocybin, we compared FC changes 1–21 days postpsilocybin to predrug baseline. The FC change analysis (described above) indicated that connectivity at the wholebrain level did not change following psilocybin (Supplementary Fig. 1). A screen was conducted with P < 0.05 threshold to identify brain networks or areas showing persistent effects. This analysis identified the anterior hippocampus as a candidate region of interest for persistent FC change (section ‘Baseline/after psilocybin FC change analysis’ in Supplementary Methods).
We assessed change in anterior hippocampus ‘FC change’ pre versus postpsilocybin using the LME model described previously. In this model, all sessions before psilocybin (irrespective or crossover order) were labelled as prepsilocybin and all sessions within 21 days after psilocybin were labelled as postpsilocybin.
As a control, we tested anterior hippocampus FC change pre versus postMTP using both the LME model, and an equivalence test. To control for potential persistent psilocybin effects, only the block of scans immediately before and after MTP were used (for example, if a participant took MTP as drug 1, then all baseline scans were labelled as ‘preMTP’ and all scans between drugs 1 and 2 were labelled ‘postMTP’).
Equivalence testing (to conclude no change in anterior hippocampus after MTP) was accomplished by setting δ = 0.5 standard deviation of FC change across preMTP sessions. We computed the 90% CI of change in FC change between pre and postMTP sessions. If the bounds of the 90% CI were within ±δ, then equivalence was determined^{123}.
Other datasets
Raw fMRI and structural data published previously^{55,56} were run through our inhouse registration and processing pipeline described above. These datasets were used for replication, external validation and generalization to another classic psychedelic (that is, LSD) for the measures described above (for example, NGSC and the MDSderived psilocybin FC dimension, dimension 1).
Using the data from ref. ^{55}: n = 15 healthy adults (five women, mean age 34.1 years, s.d. 8.2) completed two scanning sessions (psilocybin and saline) that included an eyesclosed restingstate BOLD scan for 6 min before and following i.v. infusion of drug. fMRI data were acquired using a gradientechoplanar imaging sequence, TR and TE of 3,000 and 35 ms, fieldofview 192 mm, 64 × 64 acquisition matrix, parallel acceleration factor of 2 and 90° flip angle.
Using the data from ref. ^{56}: healthy adults completed two scanning sessions (LSD and saline), which included an eyesclosed restingstate BOLD scan acquired for 22 min following i.v. drug infusion lasting 12 min. n = 20 participants completed the protocol, but data were used for n = 15 (four women; mean age 30.5, standard deviation 8.0) deemed suitable for BOLD analyses. fMRI data were acquired using a gradientechoplanar imaging sequence, TR and TE of 2,000 and 35 ms, fieldofview 220 mm, 64 × 64 acquisition matrix, parallel acceleration factor of 2, 90° flip angle and 3.4 mm isotropic voxels.
The ABCD database restingstate functional MRI^{59} (annual release v.2.0, https://doi.org/10.15154/1503209) was used to replicate the effects of stimulant use on FC. Preprocessing included framewise censoring with a criterion of frame displacement less than or equal to 0.2 mm in addition the standard predefined preprocessing procedures^{124}. Participants with fewer than 600 frames (equivalent to 8 min of data after censoring) were excluded from the analysis. Parcelwise groupaveraged FC matrices were constructed for each participant as described above for 385 regions on intertest in the brain.
Use of a stimulant (for example, MTP, amphetamine salts, lisdexamfetamine) in the last 24 h was assessed by parental report. Participants with missing data were excluded. Regression analysis was used to assess the relationship between FC (edges) and stimulant use in the last 24 h. Framewise displacement (averaged over frames remaining after censoring) was used as a covariate to account for motionrelated effects. The tvalues that reflect the relationship between stimulant use and FC were visualized on a colour scale from −5 to +5 to provide a qualitative information about effect of stimulant use on FC.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.