^a Department of Psychiatry, University of Oxford, UK
^b Centre for Mathematics of Precision Healthcare, Imperial College London, UK
^c Department of Mathematics, Imperial College London, UK
^d Centre for Psychedelic Research, Department of Brain Sciences, Imperial College London, UK
^e Department of Psychology, Yale University, USA
^f Mathematical Institute, University of Oxford, UK
^g Center for Brain and Cognition, Universitat Pompeu Fabra, Spain
^h Institució Catalana de la Recerca i Estudis Avançats (ICREA), Universitat Pompeu Fabra, Spain
ⁱ Centre for Music in the Brain, Department of Clinical Medicine, Aarhus University, Denmark
^j Institut d'Études Avancées de Paris, France
^k Life and Health Sciences Research Institute (ICVS), School of Medicine, University of Minho, Portugal
^* Corresponding author: Department of Psychiatry, University of Oxford, UK
E-mail address: joana.cabral@psych.ox.ac.uk (J. Cabral).

Table of Contents:
Abstract
1. Introduction
2. Methods
3. Participants
4. Ethics statement
5. Experimental protocol
6. Neuroimaging data acquisition & processing
 6.1. Anatomical scan acquisition
 6.2. fMRI acquisition
 6.3. fMRI processing
7. Dynamic BOLD phase-locking analysis
 7.1. Leading eigenvector of the phase-locking matrix
8. BOLD phase-locking states
 8.1. Detection of recurrent BOLD phase-locking patterns
 8.2. Probability of occurrence and switching profiles of PL states
 8.3. Selection of the optimal number of PL states
 8.4. Comparison with resting-state networks
 8.5. Global order and metastability
9. Data and code availability statement
10. Results
 10.1. Repertoire of PL states reveals canonical resting-state networks
 10.2. Consistency of between-condition differences across partition models
 10.3. Network-specific neuromodulation
 10.4. Validation using the placebo dataset
 10.5. Dynamical exploration in a low-dimensional manifold
 10.6. Trajectories between BOLD phase-locking states
 10.7. Subjective intensity of the psychedelic experience correlates with network occupancy
 10.8. Psilocybin modulates the global order of BOLD phases
 10.9. Effect of motion on PL state probabilities
11. Discussion
Acknowledgments
References

Abstract

Growing evidence from the dynamical analysis of functional neuroimaging data suggests that brain function can be understood as the exploration of a repertoire of metastable connectivity patterns ('functional brain networks'), which potentially underlie different mental processes.
The present study characterizes how the brain's dynamical exploration of resting-state networks is rapidly modulated by intravenous infusion of psilocybin, a tryptamine psychedelic found in "magic mushrooms". We employed a data-driven approach to characterize recurrent functional connectivity patterns by focusing on the leading eigenvector of BOLD phase coherence at single-TR resolution. Recurrent BOLD phase-locking patterns (PL states) were assessed and statistically compared pre- and post-infusion of psilocybin in terms of their probability of occurrence and transition profiles. Results were validated using a placebo session.
Recurrent BOLD PL states revealed high spatial overlap with canonical resting-state networks. Notably, a PL state forming a frontoparietal subsystem was strongly destabilized after psilocybin injection, with a concomitant increase in the probability of occurrence of another PL state characterized by global BOLD phase coherence. These findings provide evidence of network-specific neuromodulation by psilocybin and represent one of the first attempts at bridging molecular pharmacodynamics and whole-brain network dynamics.

1. Introduction

Brain dynamics can be understood as the exploration of activity configurations over both space and time (Tononi and Edelman, 1998; Watanabe et al., 2014; Cabral et al., 2017b; Gu et al., 2018). This exploration may be defined in terms of trajectories within a fixed repertoire of metastable activity patterns, which potentially underlie different brain processes. Indeed, a robust repertoire of large-scale functional networks has been consistently detected across individuals and neuroimaging modalities, not only during task performance but also during rest (Damoiseaux et al., 2006; De Luca et al., 2006; Mantini et al., 2007; Seeley et al., 2007; Musso et al., 2010; Brookes et al., 2011; Yeo et al., 2011; Hipp et al., 2012), indicating that brain function involves the coordinated integration of information over a repertoire of spatially distributed networks of specialized brain areas (Varela, 1979; Bressler and Menon, 2010; Menon, 2011; Cavanna et al., 2017; Lord et al., 2017). While the mechanisms driving the spontaneous formation and dissolution of functional networks remain under debate (Cabral et al., 2017a), brain function has recently been explored in terms of transitions between recurrent states of functional connectivity (Hansen et al., 2015; Cabral et al., 2017a; Cavanna et al., 2017; Tewarie et al., 2019). Although the functional relevance of these explorative dynamics during rest remains unclear (Christoff et al., 2016), recent evidence suggests that transitions between brain states are organized in a hierarchical manner (Vidaurre et al., 2017) and relate to cognitive function (Cabral et al., 2017b).

Novel insights into the neurobiological correlates of functional network states and their relationship to behavior in both health and disease can be gained by understanding how specific psychoactive compounds modulate the relative stability of functional networks over time, as well as the transitions between them. Toward this aim, the present study investigates how psilocybin (4-phosphoryloxy-N,Ndimethyltryptamine) changes the exploration of the brain's dynamical repertoire. Psilocybin is a prodrug of psilocin (4-OH-N,N-dimethyltrayptamine) - a psychoactive ingredient found in so-called "magic mushrooms" and classical tryptamine 'psychedelic'. The subjective effects of psilocybin/psilocin include broadly unconstrained perception and cognition, hyper-associative cognition and, at higher doses, a breakdown in the perception of time, space and selfhood (Griffiths et al., 2006; Carhart-Harris et al., 2014). Psilocybin provides a promising experimental framework for linking molecular pharmacodynamics to changes in the brain's dynamical repertoire because its potent psychoactive effects are selectively due to its agonist activity at the serotonin 2A (5-HT_2A) receptor (McKenna et al., 1990; Vollenweider et al., 1998; Passie et al., 2002; González-Maeso et al., 2007) - see Beliveau et al. (2017) for the highest resolution in vivo mapping of 5-HT_2A receptor densities in the human brain. Furthermore, investigating how psilocybin modulates the exploration of functional network states over time may help better understand the functional mechanisms underlying the recently demonstrated therapeutic potential of psilocybin (and other indolealkylamines) for disorders including depression, anxiety and addiction (McKenna et al., 1990; Griffiths et al., 2006; Grob et al., 2011; Meltzer et al., 2012; Johnson et al., 2014; Carhart-Harris et al., 2016).

The first fMRI investigation of psilocybin's effects on the human brain consisted of a task-free paradigm in which healthy subjects were intravenously injected with the compound inside the scanner (or a placebo, at least 1 week apart) to characterize changes in brain activity under the influence of the drug (Carhart-Harris et al., 2012). The first dynamical analysis of this dataset found increased variance of intra-network synchrony over time for a number of canonical resting-state networks under psilocybin (Carhart-Harris et al., 2014). Somewhat consistently, a subsequent analysis of this same fMRI dataset using sliding-window correlations within a specific network of four brain regions (comprising the bilateral hippocampus and anterior cingulate cortex) revealed both a larger repertoire of functional motifs under psilocybin as well as greater entropy of the motif sequence (Tagliazucchi et al., 2014). Another experiment using MEG revealed increased Lempel-Ziv complexity -a measure of neural signal diversity calculated at the single-channel level (Lempel and Ziv, 1976)- following psilocybin administration relative to placebo (Schartner et al., 2017). However, it is important to clarify that an increase in local complexity does not necessarily imply higher randomness at the larger scale. This seeming paradox can be explained by the findings in a recent study (Atasoy et al., 2018), where higher coherence - expressed as increased power and energy of connectome harmonic brain states under psilocybin - was accompanied by an enlarged repertoire of brain states, leading to higher spatial and temporal variability. Additionally, other works have shown that the serotonergic psychedelics psilocybin and LSD induce an alternative type of functional integration characterized by greater global integration (Petri et al., 2014; Roseman et al., 2014; Tagliazucchi et al., 2016). However, despite the consistent advancements, how psilocybin modulates the relative occupancy of specific functional networks over time and how it relates to the subjective psychedelic experience has not yet been explored.

Recently, a number of methodological approaches have been proposed to analyze BOLD connectivity dynamics at high temporal resolution (i.e., single volume/TR), focusing either on BOLD co-activation patterns (Tagliazucchi et al., 2012; Liu and Duyn, 2013; Karahanoglu and Van De Ville, 2015), or on BOLD phase coherence patterns (Glerean et al., 2012; Hellyer et al., 2015; Cabral et al., 2017b). While co-activation approaches (in their variant forms) are only sensitive to simultaneity in the data, phase coherence techniques can, by definition, capture temporally delayed relationships, which may be more sensitive to capture the ultra-slow oscillatory processes governing the formation of resting-state networks - as proposed by recent experimental and computational studies (Deco et al., 2009; Cabral et al., 2014a, 2017a; Ponce-Alvarez et al., 2015; Deco and Kringelbach, 2016; Gutierrez-Barragan et al., 2018; Roberts et al., 2019).

To identify recurrent patterns of BOLD phase coherence across subjects and quantify differences in the exploration of the repertoire of functional networks in the current dataset, we employed a recently developed data-driven approach, the Leading Eigenvector Dynamics Analysis (LEiDA), which captures instantaneous phase-locking patterns with reduced dimensionality by considering only the relative phase of BOLD signals (i.e., how all BOLD phases project into their leading eigenvector at each discrete time point) (Cabral et al., 2017b; Figueroa et al., 2019). LEiDA appears as a valuable step in the search for functional network recurrences in dynamical analysis because the reduced dimensionality (from an NxN matrix to a 1xN vector) allows for better convergence of the clustering algorithm, revealing robust BOLD phase-locking patterns that consistently reoccur across fMRI scans for different subjects. However, despite returning meaningful functional subsystems in previous studies (Cabral et al., 2017b; Figueroa et al., 2019), the recurrent patterns identified by LEiDA have not yet been qualitatively compared to canonical resting-state networks (Damoiseaux et al., 2006; Yeo et al., 2011).

In the current work, we therefore employed the LEiDA approach to identify recurrent BOLD phase-locking patterns (PL states) and quantified differences in their probability of occurrence and transition profiles before and after the infusion of psilocybin, while subjects were inside the MRI scanner. The validity of the results was subsequently verified using the placebo dataset, and the repertoire of PL patterns returned by LEiDA was compared to well-established resting-state networks described in the literature.

2. Methods

The full description of the participants, protocol and acquisition is provided in Carhart-Harris et al. (2012). Only a brief description is included below.

3. Participants

Employing rigorous standards for data quality, namely minimal motion during the scan, at least 21 years of age, no personal or family history of a major psychiatric disorder, no substance dependence, no cardiovascular disease, and no history of adverse response to a psychedelic drug yielded datasets for nine participants. All subjects had used psilocybin at least once before, but not within 6 weeks of the study.

4. Ethics statement

The study was approved by a National Health Service research ethics committee and all participants gave informed consent to participate in the study.

5. Experimental protocol

All subjects underwent two 12 min eyes-closed resting state fMRI scans over separate sessions, at least 7 days apart. In each session, subjects were injected intravenously with either psilocybin (2 mg dissolved in 10 ml saline, 60 s intravenous injection) or a placebo (10 ml saline, 60 s intravenous injection) in a counterbalanced design (see Fig. 2 for an illustration of the scanning paradigm). Injections were given manually by a Physician within the scanning suite. The infusions began exactly 6 min after the start of the 12 min scans and lasted 60 s. The subjective effects of psilocybin were felt almost immediately after injection and sustained for the remainder of the scanning session. This experimental approach thus provided four distinct fMRI recordings (6 min each) corresponding to each of four experimental conditions: pre/post placebo injection and pre/post psilocybin injection. To account for potential confounds related to the infusion event and to focus analyses on the active drug state, the last 5 min of each scan were used for the relevant analyses. In order to have BOLD timeseries of equal length pre- and post-injection, only the first 5 min of scanning (pre-infusion) were considered in the analysis. At the end of each session, subjects were asked to rate the overall intensity of their subjective experience under the drug (or placebo) and to comment on their wakefulness level throughout the scan. As expected, all participants rated the subjective effects of psilocybin (mean intensity = 6.9/10 ± 2.6) as much stronger than placebo (mean intensity = 0.4/10 ± 0.6) and none of the subjects reported falling asleep during either scanning session.

6. Neuroimaging data acquisition & processing

6.1. Anatomical scan acquisition

Neuroimaging data were acquired using a 3T GE HDx MRI system. Anatomical scans were performed before each functional scan and thus prior to administering either the drug or placebo. Structural scans were collected using a 3D fast spoiled gradient echo scans in an axial orientation, with field of view = 256 × 256 × 192 and matrix = 256 × 256 × 192 to yield 1 mm isotropic voxel resolution (repetition time/echo time TR/TE = 7.9/3.0 ms; inversion time = 450 ms; flip angle = 20).

6.2. fMRI acquisition

BOLD-weighted fMRI data were acquired using a gradient echo planar imaging sequence, TR/TE 3000/35 ms, field-of-view = 192 mm, 64 × 64 acquisition matrix, parallel acceleration factor = 2, 90 flip angle. Fifty-three oblique axial slices were acquired in an interleaved fashion, each 3 mm thick with zero slice gap (3 × 3 × 3 mm voxels).

6.3. fMRI processing

fMRI data were processed using MELODIC (Multivariate Exploratory Linear Decomposition into Independent Components) (Beckmann and Smith, 2004), part of FSL (FMRIB's Software Library, www.fmrib.ox.ac.uk/fsl). The default parameters of this imaging pre-processing pipeline were used on all participants: motion correction using MCFLIRT (Jenkinson et al., 2002), non-brain removal using BET (Smith, 2002), spatial smoothing using a Gaussian kernel of FWHM 5 mm, grand-mean intensity normalization of the entire 4D dataset by a single multiplicative factor and linear de-trending over 50 s intervals.

We used the Anatomical Automatic Labeling (AAL) atlas (Tzourio-Mazoyer et al., 2002) to parcellate the MNI brain into N = 90 cortical and sub-cortical non-cerebellar brain areas and the BOLD signals were then averaged over all voxels belonging to each brain area using FSL tools. The BOLD signals in each of the 90 brain areas were subsequently band-pass filtered between 0.02 and 0.1 Hz (using a 2nd order Butterworth filter), discarding in this way the high frequency components associated to cardiac and respiratory signals (>0.1 Hz), and focusing on the most meaningful frequency range of resting-state fluctuations (Biswal et al., 1995).

For each subject, this procedure was applied separately to the fMRI data from all four experimental conditions listed above, resulting in four NxT BOLD datasets where N = 90 is the number of brain areas and T = 100 is the number of TRs.

7. Dynamic BOLD phase-locking analysis

To compute the phase alignment between each pair of AAL regions, first the BOLD phases, θ(n,t), were estimated using the Hilbert transform for each BOLD regional timecourse (see Fig. 1A) (Glerean et al., 2012; Cabral et al., 2017b). The Hilbert transform expresses a given signal x as x(t) = A(t)*cos(θ(t)), where A is the time-varying amplitude and θ is the time-varying phase. In Fig. 1A, we represent the BOLD signal phase of one area n over time as e^iθ(t) with sin(θ(t)) representing the imaginary part of the analytic phase, and cos(θ(t)) representing its real part (black dotted lines). We show that cos(θ(t)) captures the oscillatory dynamics of the original BOLD signal (green) - which is mostly preserved after band-pass filtering between 0.02 and 0.1 Hz (blue) - but with constant amplitude (i.e., between -1 and 1). The arrows in red represent the Hilbert phases at each TR, which can be projected into the complex plane (i.e., the plane defined by the real and imaginary axes, represented at t = 0).

Fig. 1B shows, at a single time point t, all N = 90 BOLD phases in the cortical space, placed at the center of gravity of each brain area. Here, the arrows are colored according to their direction when projected into the leading eigenvector of phase coherence. On the middle of panel B, the same N = 90 BOLD phases are plotted in the complex plane, i.e., all centered at the same origin.

To obtain a whole-brain pattern of BOLD phase coherence at each single time point t, we compute a dynamic Phase-Locking matrix dPL(n,p,t) which estimates the phase alignment between each pair of brain areas n and p at each time t using Equation (1):

dPL (n, p, t) = cos(θ(n, t) - θ(p, t)).            (1)

Using the cosine function, two areas n and p with temporarily aligned BOLD signals (i.e., with no phase difference) at a given TR will have a phase-locking value dPL(n, p, t) = cos(0°) = 1. On the other hand, time points when the BOLD signals have180° phase difference (in complex plane) will have dPL(n, p, t) = cos(180°) = -1. The cosine function being even, the dPL(t) matrix is symmetric with values ranging between -1 and 1.

The resulting dPL for each subject in each condition is thus a three-dimensional tensor with size NxNxT, where N = 90 is the number of brain areas and T = 100 is the total number of time points.

Fig. 1. Detection of recurrent BOLD phase-locking (PL) patterns. A) The BOLD signal in a given brain area n (green) is first band-pass filtered between 0.02 and 0.1 Hz (blue) and then transformed into an analytic signal (with real and imaginary components) using the Hilbert transform. The phase dynamics of the analytic signal can be represented over time by eiθ (black line), where the real part is captured by cos(θ) and the imaginary part by sin(θ) (black dotted lines). The red arrows represent the BOLD phases at each TR. B) At a single time point, the BOLD phases in all N = 90 brain areas are represented both in the cortical space (arrows centered at the center of gravity of each brain area) and in the complex plane (i.e., the unit circle with real and imaginary axes, where all phases are centered at the same origin). The phase-locking matrix at time t captures the BOLD phase alignment between each pair of brain areas. (Right) The leading eigenvector of the phase-locking matrix at time t, V1(t), is the vector that best captures the main orientation of all BOLD phases (represented as a black dashed arrow in the unit circle in the left), where each element in V1(t) corresponds to the projection of the BOLD phase in each area into V1(t). The elements of V1(t) are colored according to their relative direction with respect to V1 (red: positive; blue: negative). C) We take the leading eigenvector V1(t) as a low-dimensional representation of the BOLD phase-locking patterns over time. To identify recurrent phase-locking patterns, we apply a clustering algorithm (k-means), which divides the pool of data points into a predefined number of clusters k. Each cluster is represented by a central vector, which we take to represent a recurrent BOLD phase-locking pattern, or PL state.

7.1. Leading eigenvector of the phase-locking matrix

To characterize the evolution of the dPL matrix over time with reduced dimensionality, we employed a method termed Leading Eigenvector Dynamics Analysis (LEiDA) (Cabral et al., 2017b; Figueroa et al., 2019). The leading eigenvector of the NxN phase-locking matrix at time t, V₁(t), is a Nx1 vector that captures the main orientation of BOLD phases over all areas, where each element in V₁(t) represents the projection of the BOLD phase in each brain area into the leading eigenvector (Fig. 1B, right). When all elements of V₁(t) have the same sign, all BOLD phases are pointing in the same direction (half-plane) with respect to the orientation determined by V₁(t), which is indicative of a global mode governing all BOLD signals. If instead the first eigenvector V₁(t) has elements of different signs (i.e., positive and negative), the BOLD signals follow different directions with respect to the leading eigenvector, which naturally divides the brain areas into 2 clusters according to their BOLD phase relationship (see Fig. 1B). Moreover, the magnitude of each element in V₁(t) indicates the 'strength' with which brain areas belong to the communities in which they are placed (see Newman (2006) and the Supplementary Information for further details). Since V and -V span the same one-dimensional subspace, we use a convention ensuring that most of the elements have negative values, because, as we will show, our results revealed that the smallest group of areas whose BOLD phases do not follow the global mode reveal meaningful functional brain networks (i.e., canonical resting-state networks).

This approach substantially reduces the dimensionality of the data, while still explaining most of the variance of BOLD phase coherence (see the Supplementary Fig. S1 where we show that the leading eigenvector consistently represents >50% of the variance in phase coherence at all time points).

8. BOLD phase-locking states

8.1. Detection of recurrent BOLD phase-locking patterns

To identify recurrent PL patterns, we applied a k-means clustering algorithm to divide the set of 1800 eigenvectors (corresponding to all 100 TRs of all 9 subjects both before and after psilocybin injection) into a predefined number of clusters k (see Fig. 1C), with higher k revealing more rare and more fine-grained patterns. Since the optimal number of functional networks to consider remains an open question, we ran the k-means clustering algorithm with k ranging from 5 to 10 (i.e., dividing the set of eigenvectors into k = 5, 6, ..., 10 clusters) to cover the range of functional networks commonly reported in the resting-state fMRI literature (Beckmann et al., 2005; Damoiseaux et al., 2006; Yeo et al., 2011).

For each partition model considered (i.e., k = 5 to k = 10), the clustering returns k cluster centroids in the shape of Nx1 vectors V_C, which represent the average vector of each cluster (in Fig. 2A we show the 7 central vectors obtained with k = 7; the centroids obtained for all solutions ranging from k = 5 to k = 10 are shown in Supplementary Fig. S2). We take these central vectors as representing recurrent BOLD phase-locking states, or PL states. As shown in Fig. 2B, each PL state can be represented as a network in cortical space, where the value of V_C(n) is used to scale the color of each brain area and links are plotted between areas with positive sign to highlight the network detaching from the global mode (in Fig. 2 we used V_C(n) > 0.3 for visualization purposes only). Also, to facilitate visualization and interpretation of PL states, the cluster centroid vectors V_C can be rendered onto a cortical surface, e.g. using the HCP Workbench (Fig. 2C). Finally, we note that the PL states can also be represented back into matrix format (NxN) by computing the outer product V_C. V_C^T resulting in a matrix of rank 1 (i.e., a matrix that is obtained from a single vector), with positive values between all elements with the same sign in V_C (be they positive or negative), and negative values between elements of different signs (Fig. 2D).

8.2. Probability of occurrence and switching profiles of PL states

The clustering assigns to each TR a single PL state, by selecting the closest centroid V_C at each TR, as illustrated by the color-coded bars in Fig. 3A for one representative subject (the color-coded cluster time courses are shown for all 9 subjects in Supplementary Figs. S3 and S4). Using the state time courses, we calculated the probability of occurrence of each state, which is simply the number of epochs assigned to a given PL state divided by the total number of epochs (TRs) in each scanning session (see Fig. 3A for an illustration with k = 7). For each partition model (i.e., with k = 5 to k = 10) the probabilities of each PL state were calculated for each subject before and after psilocybin injection separately.

Differences in probabilities of occurrence before and after injection were statistically assessed using a permutation-based paired t-test. This non-parametric test uses permutations of group labels to estimate the null distribution, which is computed independently for each experimental condition (before versus after psilocybin injection). For each of 1000 permutations, a t-test is applied to compare populations and a p-value is returned.

Fig. 2. Repertoire of recurrent BOLD phase-locking patterns (PL states) obtained by clustering the leading eigenvectors of BOLD phase coherence (here with k = 7). Seven recurrent PL states were obtained from unsupervised clustering of the 1800 eigenvectors of BOLD phase coherence over time, sorted (left to right) according to decreasing probability of occurrence before psilocybin injection. Each PL state is represented by the central vector of each cluster VC as: A) a bar plot showing the N elements in VC, representing the projection of the BOLD phase in each brain area into the leading eigenvector; B) a network in cortical space, where the value of VC(n) is used to scale the color of each brain area and links are plotted between areas with VC(n) > 0.3 to highlight the network detaching from the global mode; C) cortical renderings in AAL space of regional contributions to the respective VC and D) a matrix obtained by calculating the outer product of VC, where positive values are the product of VC elements with the same sign, be they positive or negative.

8.3. Selection of the optimal number of PL states

In the current study we did not aim to determine the optimal number of PL states governing resting-state activity, but rather to investigate if there are PL state(s) that significantly differ in their probability of occurrence following the administration of psilocybin. As such, we employed a data mining approach to detect meaningful activity patterns (PL states) from the session during which psilocybin was infused, and left the placebo session (recorded more than 7 days apart) as a validation dataset (see section Validation using the Placebo dataset).

To search for PL patterns affected by psilocybin, we varied k (number of clusters) between 5 and 10, and for each k examined how the probability of occurrence of each PL state changed after the injection of psilocybin. We subsequently analyzed the robustness of the results over the range of partition models by plotting the p-values with respect to different significance thresholds, i.e., the standard α₁ = 0.05 or the Bonferroni corrected significance threshold α₂ = 0.05/k to correct for the number of independent hypotheses compared in each partition model.

Fig. 3. Assignment of one BOLD Phase-Locking (PL) state to each TR and corresponding probabilities of occurrence for each condition and participant. A) For each of the 9 participants, 2 fMRI sessions were recorded at least one week apart. In the middle of each session, either Psilocybin or a Placebo was intravenously administered inside the scanner. The color-coded bars behind the BOLD signals represent the PL state assigned to each TR (corresponding to the 7 PL states reported in Fig. 2). Results are shown here for subject 9 (the remaining participants are shown in SI Figs. S3 and S4). B) The probabilities of occurrence of each PL state were obtained for each of the 4 conditions and for each subject separately, which were used for the subsequent statistical analysis.

For the subsequent validation using the placebo condition, we used exactly the same PL states obtained previously for the psilocybin session for the selected optimal k, and evaluated how these same PL states were explored during the placebo session. As in the previous analysis, we first obtained the 1800 leading eigenvectors V₁ (9 subjects, 2 conditions, 100 TRs each) and ran a single iteration of the k-means algorithm with k = 7, defining as 'start vectors' the 7 cluster centroids V_C obtained before. This computes the cosine distance between each observation in the placebo dataset and the previously defined 7 cluster centroids, returning a sequence of cluster time courses for the placebo fMRI session. Critically, this approach allows a direct comparison with the psilocybin fMRI sessions recorded several days apart, enabling us to verify whether the probabilities of approaching a given cluster centroid in the placebo sessions were: i) consistent with the pre-psilocybin condition and ii) different from the post-psilocybin condition.

In a follow-up analysis, we compared the number of within-subject occurrences of the most meaningful PL states before and after the psilocybin injection using a paired t-test.

Fig. 4. Overlap of the 6 non-global PL states with seven canonical resting-state networks (RNSs). A) The 7 cortical networks defined in Yeo et al. (2011) were transformed into AAL space by counting the number of 2 mm3 MNI voxels in each AAL brain area that belong to each of the 7 networks, obtaining in this way a vector with N = 90 elements that we can compare with the PL states obtained in this work. B) Representation in cortical space of the vectors shown in (A) colouring in red the spheres representing the brain areas contributing to each RSN. Red links were plotted between AAL areas with >400 MNI voxels contributing to each RSN. C) We compared the 6 non-global PL states obtained with LEiDA (for k = 7) with the 7 RSNs by computing the Pearson's correlation between the 1 × 90 centroid vectors VC (shown in Fig. 3) and the vectors shown in (A). Asterisks denote significant correlations with p-value < 0.001.

8.4. Comparison with resting-state networks

We first transformed the 7 RSNs defined in 2 mm³ MNI space by Yeo et al. (2011) into 7 vectors with 90 elements each (shown in Fig. 4), where each element scales the contribution of each AAL brain area to the corresponding RSN. Subsequently, we computed the bivariate correlation with the centroid vectors V_C, setting all negative elements in V_C to zero, keeping only the values of the positive elements (the red links in Fig. 2B).

8.5. Global order and metastability

Using a perspective from dynamical systems theory, we investigated how psilocybin affected the order and stability of BOLD signals. At each instant of time, the degree of order between BOLD phases θ(n, t), n = 1,..,90, can be quantified using the magnitude of the Kuramoto Order Parameter, OP(t), which can range between 0 (when all BOLD signals are out of phase) and 1 (when all BOLD signals are in phase):

$$\textit{OP(t)} = \frac{1}{N}\vert\sum_{n=1}^N e^{iθ(n,t)}\vert\text{.}$$

The mean magnitude of the Order Parameter over time informs whether the system is mostly incoherent, partially synchronized or fully synchronized (Acebrón et al., 2001; Acebrón et al., 2005; Cabral et al., 2011; Deco and Kringelbach, 2016). Moreover, the standard deviation of the order parameter, STD(OP), can be used to characterize the degree of metastability in the system (Shanahan, 2010; Cabral et al., 2014b). In other words if the order parameter is constant over time, it means that the system is in a stable equilibrium, be it synchronized or not. When more than one weakly stable states co-exist and the system systematically switches from one to another, then the variance of the order parameter increases, indicating a dynamical system with metastability. To investigate the effects of psilocybin on the stability of the global order of BOLD phases, we therefore calculated the standard deviation of the order parameter over time within each subject before and after the psilocybin infusion. Between-condition differences on these measures were assessed using a one-sample t-test.

9. Data and code availability statement

fMRI data from the 9 participants in the 4 conditions is available on request. Please contact Robin L Carhart-Harris r.carhart-harris@imperial.ac.uk at the Centre for Psychedelic Research, Department of Brain Sciences, Imperial College London.

The Matlab codes developed for the analysis are made available open source at
github.com/juanitacabral/LEiDA_Psilocybin, together with the post-processed data.

10. Results

10.1. Repertoire of PL states reveals canonical resting-state networks

The repertoire of PL states obtained when dividing the set of leading eigenvectors into seven clusters is plotted in Fig. 2, where each state is represented by its central vector V_C as: (A) a bar plot showing the projection of the BOLD phase in each brain area into the leading eigenvector, (B) a network in cortical space, where links are plotted between areas with V_C(n) > 0.3 to highlight the network detaching from the global mode; (C) a rendering of regional values of V_C onto the cortical surface, and (D) a matrix obtained by calculating the outer product of V_C. While partitions with different k were produced, we selected the partition into k = 7 PL states because it returned the PL states that most significantly differed in terms of probability of occurrence after psilocybin injection, as we will explain in section Results - Consistency of between-condition differences across partition models.

Results reveal that the most probable BOLD phase-locking state (PL state 1) corresponds to a state where all BOLD signals are following one main direction (i.e., all projecting toward the same direction into the leading eigenvector), in line with findings from prior studies using LEiDA (Cabral et al., 2017b; Figueroa et al., 2019). While this global modulation occurs for a significant fraction of the time, we detect a number of non-global patterns that transiently and recurrently shape the phase alignment of BOLD signals during the remaining time.

The non-global PL states 2-7 shown in Fig. 3 reveal striking spatial similarities with canonical resting-state networks reported in the literature. We verified this overlap by computing the correlation between the 7 resting-state networks (RSNs) defined in Yeo et al. (2011) and each PL state obtained herein (Fig. 4). Corroborating our observations, we found that all non-global PL states show strong and statistically significant spatial overlap with well-documented resting-state networks, namely PL state 6 reveals a strong correlation with the Visual network (Pearson's r = 0.79, p = 2.3*10^-20), PL state 3 correlates with the Frontoparietal network (r = 0.70, p = 1.1*10^-14), PL state 5 correlates with the Somatomotor network (r = 0.65, p = 6.1*10^-12), PL state 4 correlates with the Limbic network (r = 0.36, p = 5.5*10^-4) and PL state 2 correlates with the Default Mode network (r = 0.35, p = 5.5*10^-4). Finally, PL state 7 reveals approximately equal contributions of the Ventral Attention (r = 0.43, p = 2.0*10^-5) and the Somatomotor networks (r = 0.43, p = 2.8*10^-5).

10.2. Consistency of between-condition differences across partition models

The repertoire of PL states obtained depends on the number of clusters k defined in the k-means clustering algorithm, with a higher k generally revealing more fine-grained, less frequent and often less symmetric networks. Here, we chose to cover a range of partition models (i.e., with k varying gradually between k = 5 to k = 10) and examined the PL configurations that most significantly differentiate brain activity under psilocybin in each partition model relative to the pre-injection baseline.

In Fig. 5, we show, for each partition into k PL states (horizontal axis), the p-values obtained from between-condition comparisons (pre vs post-psilocybin) in terms of probability of occurrence of each of the corresponding k PL states. Since a higher number of comparisons increases the chances of false positives (i.e., of p-values falling below the standard threshold of α₁ = 0.05, shown in red in Fig. 5C), we Bonferroni correct the threshold by the number of independent hypotheses tested within each partition model, i.e., α₂ = 0.05/k (green dashed line in Fig. 5C). Note that across partition models, the null hypotheses being tested are not independent from each other (as shown in the correlations reported in Fig. 5D), so correcting for the full number of comparisons would not be adequate here.

We find that, for the entire range of partition models explored (i.e., 5 ≤ k ≤ 10), the PL state that most strongly alters its probability of occurrence after the psilocybin infusion is a frontoparietal network (Fig. 5A) passing the corrected threshold α₂ for all k except for k = 6 (p = 0.0095), which despite falling below α₁, narrowly misses the α₂ cutoff. The most significantly different PL state of each partition model is reported in Fig. 5B in vector format. These PL states were highly correlated across clustering solutions (CC > 0.86 for all pairs of PL states, Fig. 5D), which indicates that they refer to the same underlying functional network, with small differences arising from the number of output states constrained by k. Furthermore, we note that a second network exhibited significant differences between conditions surviving the correction for the number of clusters (α₂ = 0.05/k) for the clustering solutions: k = 5 (p = 0.0078) and k = 7 (p = 0.0068) (Fig. 5C, second smallest p-value for each k). This network corresponds to the globally coherent PL state 1, which also showed a high level of consistency across the range of clustering solutions, with some brain areas being more strongly aligned to the global mode of BOLD phase coherence than others (see Supplementary Fig. S2).

Fig. 5. Detection of the PL states that most significantly change their probabilities of occurrence after psilocybin infusion across partition models. A) The PL state differing the most between conditions is a functional network consisting of lateral frontoparietal areas, shown here for k = 7 (PL state 3 from Fig. 2). B) Bar plots showing, for each partition into k = 5 up to k = 10 PL states, the cluster centroid of the PL state that most significantly changed in probability of occurrence between conditions (p-values reported in the title are colored according to significance level). C) For each partition model, we plot the p-values associated to the comparison between all PL state probabilities before and after psilocybin injection. We find that, although most PL states do not show significant differences between groups (black dots falling above the α1 = 0.05 threshold), some PL states show significant differences in probability of occurrence between groups, falling below (or very near for k = 6) the corrected threshold by the number of clusters (α2 = 0.05/k, green dashed line). D) The PL states representing the Frontoparietal network in different partition models (shown in B) are highly correlated with one another.

10.3. Network-specific neuromodulation

The probabilities of occurrence of each PL state in each of the 4 experimental conditions for the optimal solution with k = 7 are shown in the bar plots of Fig. 6B. As predicted from the previous analysis, we find that for k = 7, psilocybin significantly alters the occurrence of two of the seven BOLD phase-locking modes. The probability of occurrence of PL state 3 (red) significantly decreased from 14.3 ± 2.4% pre-injection (third most visited state) to 4.1 ± 1.2% following the psilocybin administration (p = 0.00034 uncorrected, p' = 0.0021 after correction). Conversely, a significant increase in the probability of occurrence of PL state 1 (dark blue) was also observed, shifting from 26.7 ± 4.5% pre-psilocybin injection to 37.2 ± 6.0% post-psilocybin injection (p = 0.0068 uncorrected, p' = 0.047 after correction). In contrast, the probability of occurrence of all other PL states remained statistically similar before and after the psilocybin infusion (all p-values > 0.05).

Fig. 6. Psilocybin modulates the exploration of BOLD Phase-Locking configurations. A) Repertoire of PL states detected with LEiDA for k = 7 from the session when psilocybin was infused (same as shown in Fig. 2). B) Probability of occurrence (mean ± standard error of the mean across subjects) of each PL state in each of the four experimental conditions (pre/post psilocybin and pre/post placebo injection). We find that after the infusion of psilocybin, PL state 3 (red) significantly decreases in probability of occurrence from an average of 14.3-4.1% (p = 0.00034), whereas PL state 1 (blue) significantly increases after psilocybin injection (p = 0.0068). All other PL states did not change in fractional occupancy with psilocybin. Importantly, when the BOLD phase leading eigenvectors from the placebo session (900 before and 900 after infusion) were matched to the 7 PL cluster centroids obtained from the psilocybin session, we found that all the probabilities of occurrence were consistent with the resting state before psilocybin injection (despite being more than a week apart), and exhibited the same significant differences with respect to the psilocybin condition. Asterisks denote the presence of statistically significant differences (* < 0.05; ** < 0.05/k correcting for the number of independent hypotheses tested). C) Representation of all PL configurations in each experimental condition on a low-dimensional manifold highlighting the reduced number of occurrences of the Frontoparietal network (red dots) after psilocybin injection. In each scatter plot, each dot represents one PL configuration, captured by V1(t). Each panel corresponds to a different session (before/after psilocybin/placebo). Dots are colored according to the closest cluster centroid vector, VC, using the same color code as in B. It can clearly be observed that the number of red dots is substantially reduced after the psilocybin injection. This representation serves to illustrate how the reduced occurrence of the Frontoparietal network under psilocybin can be interpreted in terms of fewer excursions into a specific region of a low-dimensional manifold of network configurations.

10.4. Validation using the placebo dataset

The 7 PL states shown in Fig. 2 were defined solely from half of the fMRI scans in the dataset, i.e., the scans from the 9 participants recorded during the day in which psilocybin was infused (including the pre-injection baseline), consisting of 9 2 100 TRs. Using a data mining approach, we searched for the partition model that returned the most significant differences between pre- and post-injection conditions.

To further validate our findings, we verified whether these same 7 PL states were expressed in the fMRI session during which placebo was infused - recorded more than a week apart -, and if they occurred with similar probabilities as during the baseline condition before psilocybin was injected. To do so, we took exactly the same cluster centroids (shown in Fig. 2A) and evaluated their probabilities of occurrence both before and after the placebo infusion.

Validating our results, we found that the probabilities of occurrence of all 7 PL states both before and after placebo injection did not differ from the baseline condition before psilocybin was infused (all p-values > 0.05). Notably, this indicates that the dynamical exploration of PL states assessed with LEiDA remained stable across subjects under baseline conditions over a period of more than a week (Fig. 6B). Corroborating our previous findings, the fractional occupancy of PL state 3 remained within normal levels both before (15.0 ± 3.3%) and after (14.3 ± 2.9%) placebo injection (p = 0.49 uncorrected), and was in each case significantly higher than the probability of occurrence found during the psychedelic state (4.1 ± 1.2%, p = 0.00057 and p = 0.00014, respectively, both falling below the corrected threshold α₂ = 0.05/k = 0.007). Moreover, the probability of occurrence of PL state 1, remained unchanged following the placebo injection (23.6 ± 5.6% versus 27.4 ± 3.1%, p = 0.38), being substantially lower than the probability of occurrence of PL state 1 under psilocybin (37.2 ± 6.0%, p = 0.00051 before, p = 0.012 after placebo injection), with the pre-placebo comparison surviving correction for multiple comparisons with α₂.

10.5. Dynamical exploration in a low-dimensional manifold

The non-stationary dynamics of the PL leading eigenvectors over time can be interpreted as the constant exploration of a cloud of PL configurations in a multi-dimensional space (with N = 90 dimensions in the selected parcellation). In order to visualize this cloud of observations in a low-dimensional manifold (Fig. 6C), we applied a Principal Components Analysis (PCA) to all leading eigenvectors across conditions and subjects (90 × 3600) and used the first 3 Principal Components (3 × 3600) to project each observation in a low-dimensional manifold. In this way, each 90 dimensional leading eigenvector V₁ obtained at each TR can be plotted as a dot in a 3-dimensional (3-D) space defined by the first 3 principal components and colored according to its closest centroid vector, V_C, using the same color scheme from Fig. 6B. By construction, the k-means clustering algorithm produces compact clusters (i.e., non-overlapping and well separated in the 90 dimensional space where the leading eigenvectors are defined). As can be seen in Fig. 6C, retaining the first 3 principal components is sufficient to conserve well-defined clusters (see also Supplementary Fig. S3 to view the 3-D manifolds from a different angle).

This representation serves to illustrate how the occurrence of the different functional networks can be interpreted mechanistically in terms of excursions into distinct regions of a low-dimensional manifold of network configurations. The visual rendition of the data from all subjects in the 4 different experimental conditions shown in Fig. 6C highlights the reduced number of excursions into the region of the manifold associated to the frontoparietal PL state under psilocybin (number of red dots in the upper right-hand plot), compared to the 3 non-psychedelic conditions (before psilocybin and before/after placebo injection). Although less evident, it can also be seen that the blue cluster, corresponding to the globally coherent PL state 1, becomes denser (more dots) after the psilocybin injection. The statistical analysis shown in Fig. 5C serves to define the partition model that optimally delineates the region of the cloud affected after psilocybin injection.

10.6. Trajectories between BOLD phase-locking states

In order to explore the trajectories between the different PL states, we report in Fig. 7 the switching matrices containing the mean probabilities of, being in a given PL state (rows), transitioning to any of the other PL states (columns) both before (left) and after the psilocybin injection (right). Switching probabilities were calculated for each subject and each condition and then statistically compared using a permutation-based paired t-test (5000 permutations). To facilitate the interpretation of the differences found in PL state switching profiles, we provide an illustration of the transitions that were most significantly affected by psilocybin (solid red: increased; dashed blue: decreased). We find that the probability of transitioning from any given PL state to the frontoparietal PL state 3 is consistently reduced under psilocybin effects. Conversely, all PL states except one were more likely to transition toward the PL state of global coherence following the psilocybin injection. The p-values for each state-to-state transition in the switching matrix are provided in Supplementary Table ST1.

10.7. Subjective intensity of the psychedelic experience correlates with network occupancy

After each scanning session participants were asked to rate the subjective intensity of the psychedelic (or placebo) experience on a scale of 1-10. We correlated those subjective intensity ratings with the probability of occurrence of each of the 7 PL states (Fig. 8). We found that the occurrence of frontoparietal network (PL state 3) was negatively correlated with the subjective rating of the psychedelic experience (Pearson's r = -0.56; p = 0.0083, one-tailed). Conversely, a near-significant positive correlation was found between the fractional occupancy of the globally coherent PL state 1 and the subjective intensity of the psychedelic experience (r = 0.36, p = 0.07, one-tailed). Other PL states whose fractional occupancy was not changed by psilocybin did not correlate with this behavioral measure.

10.8. Psilocybin modulates the global order of BOLD phases

We also investigated how psilocybin affected the global order of BOLD phases using the Kuramoto Order Parameter (OP) and, although the mean synchrony degree did not reveal significant differences between conditions (p > 0.05), we found that the degree of whole-brain metastability (i.e., how much the synchrony degree fluctuates over time) significantly increased following the psilocybin injection compared to baseline (standard deviation, STD(OP) = 0.18 ± 0.01 vs STD(OP) = 0.16 ± 0.015, t(8) = 2.41, p = 0.04 using a paired sample t-test). In Supplementary Fig. S4, we report the OP over time before and after the psilocybin injection in all subjects.

In addition, we investigated how the OP relates with the different PL states. As expected, we found that the probability of occurrence of the PL state of global coherence (state 1) at the single subject level during a given scan is strongly correlated with the mean of the global OP (cc = 0.825, p = 2.9*10^-5) (Supplementary Fig. S4), indicating that this PL state is associated with an increase in global order, whereas all other PL states are consistently found for epochs of lower OP. In Supplementary Fig. S4, we report the global OP obtained while each of the 7 PL states was active, according to the color-code used in Fig. 6B. This allows a visual representation of all PL state assignments in all subjects at the single TR resolution, showing that PL state 1 is the most frequently visited state - both before and after the psilocybin injection - and also occurs in epochs of high OP, indicating a high level of phase coherence across the network nodes.

Finally, we ran a paired-sample t-test to confirm that, at the individual subject level, the mean number of occurrences of PL state 3 (overlapping with the Frontoparietal network) was significantly decreased following the psilocybin injection relative to the pre-injection baseline from 14.3 ± 2.5% to 4.1 ± 1.2% (t(8) = 3.88, p = 0.005). Conversely, a paired-sample t-test showed a significant increase in the mean number of occurrences of PL state 1 (globally coherent state) following the psilocybin injection relative to the pre-injection baseline from 26.% ± 4.5% to 37.2 ± 6.0% (t(8) = 2.75, p = 0.01). These statistical results reinforce the validity of the findings reported in Fig. 6B.

Fig. 8. Subjective intensity rating of the psychedelic experience is related to dynamic exploration of PL states. After the infusion of either Psilocybin (red) or Placebo (blue), participants were asked to rate the subjective intensity of their experience which in each case was far stronger in the psilocybin condition (mean intensity = 6.9 ± 2.6) than placebo (mean intensity = 0.4 ± 0.6) on a scale from 0 to 10. We further investigated the relationship between the subjective ratings and the probability of occurrence of each of the 7 PL states. We found that the occurrence of the frontoparietal network (PL state 3) was negatively correlated with the ratings on psychedelic experience (correlation coefficient cc = -0.56; p = 0.0083, one-tailed), with higher intensity ratings being associated with a lower probability of occurrence of PL state 3 (black lines link the same participant in the placebo (blue) and psilocybin (red) conditions). Conversely, a near-significant positive correlation was found between the fractional occupancy of the globally coherent PL state 1 and the subjective intensity of the psychedelic experience (cc = 0.36, p = 0.073, one-tailed).

10.9. Effect of motion on PL state probabilities

We performed follow-up analyses to confirm that the observed effects of psilocybin on the fractional occupancy of PL state 1 and PL state 3 were not due to differences in participant motion during the scan. For the nine subjects included in the study, we assessed motion by calculating the mean framewise displacement (FD) for each subject. This variable measures movement of any given frame relative to the previous frame. We then conducted two separate analyses to examine the potential effects of motion. First, a paired t-test showed that the mean FD was not significantly different between the post-placebo and post-psilocybin injection conditions: t(8) = 1.28, p = 0.24. Furthermore, the fractional occupancy of PL state 3 post-psilocybin did not correlate with the mean FD (correlation coefficient, CC = -0.14, p = 0.72) and neither did the fractional occupancy of PL state 1 post-psilocybin (CC = 0.04, p = 0.91).

11. Discussion

The present study is one of the first to investigate how the exploration of the brain's repertoire of functional networks is rapidly modulated by a psychoactive molecule, here the serotonergic psychedelic psilocybin. Specifically, we found that the dynamical exploration of a fixed repertoire of discrete functional networks at rest, defined as recurrent BOLD phase-locking patterns over time (PL states), is significantly altered after a psilocybin infusion. The strongly reduced expression of a BOLD PL state overlapping with the previously described frontoparietal control system (Vincent et al., 2008; Yeo et al., 2011) suggests that this particular functional subsystem is related to the psychoactive effects of psilocybin. Moreover, we found that the decreased expression of the Frontoparietal network under the influence of the drug only increased the stability of a state of global coherence (PL state 1) and did not change the expression of the remaining PL patterns. These new findings represent one of the first attempts at bridging the gap between molecular pharmacodynamics and the brain's network dynamics.

The highly significant disruption of a BOLD PL pattern closely overlapping with the Frontoparietal control system is consistent with prior human neuroimaging studies of psilocybin's effects. Indeed, a previous analysis of the same dataset revealed that psilocybin modified the BOLD spectral content in a distributed Frontoparietal network. Specifically, the Frontoparietal system of interest showed significantly decreased power at low frequencies (and decreased power spectrum scaling exponent) under psilocybin (Tagliazucchi et al., 2014). Similarly, an MEG study reported decreased oscillatory power under psilocybin in a bilateral Frontoparietal network derived from beta-band activity (Muthukumaraswamy et al., 2013). Findings from these prior studies indicate a rapid desynchronization of the Frontoparietal control system under psilocybin and are thus consistent with the present results. Furthermore, the disengagement of the Frontoparietal system can be interpreted in light of the cognitive and behavioral effects of psilocybin reported by the study participants. Evidence suggests the Frontoparietal system plays a key role in cognitive control by focusing attention on goal-relevant information through top-down mechanisms (Cole and Schneider, 2007; Dosenbach et al., 2008; Xin and Lei, 2015). Such goal-directed focus does indeed appear to be impaired under moderate-high doses of psilocybin (Bayne and Carter, 2018) - which would be a reasonable way of describing the present dosage. Declines in frontoparietal FC have been associated with difficulty suppressing goal-irrelevant information, and greater distractibility (Smallwood et al., 2012). These observations are broadly consistent with the reports of increased freely wandering thoughts and hyper-associative cognition provided by participants in the present cohort (Carhart-Harris et al., 2012; Turton et al., 2014).

It is interesting to consider how the known pharmacodynamics properties of psilocybin could mechanistically explain its destabilizing effect on the Frontoparietal PL state dynamics described above. The psychoactive effects of psilocybin are primarily due to its agonist activity at the serotonin 2A (5-HT_2A) receptor (Vollenweider et al., 1998; Kometer et al., 2013). The distribution and density of 5-HT_2A receptors in the human brain (and other 5-HT receptor subtypes) were recently estimated in a high-resolution atlas using molecular and structural neuroimaging data in a large cohort of healthy individuals in combination with postmortem human brain autoradiography (Beliveau et al., 2017); 5-HT_2A receptors are widely distributed across the brain but most densely expressed in the cortex. Interestingly, the 5-HT_2A receptor density map obtained by Beliveau and colleagues shows a strong qualitative overlap between 5-HT_2A receptor expression on the lateral surface of the brain, and the regions implicated in the Frontoparietal PL state detected in the present study. Moreover, an fMRI study of the synthetic 5-HT_2A receptor agonist and classical psychedelic lysergic acid diethylamide (LSD) has reported significantly increased functional connectivity density (FCD, quantified as the average functional connectivity to the rest of the brain) of a similar Frontoparietal network which, in turn, strongly overlapped with 5-HT_2A receptor concentrations (Tagliazucchi et al., 2016). The increase in FCD of the Frontoparietal network can be interpreted as reduced within-network integrity and is thus consistent with the present findings, where the PL state 3 (high within-network coherence and low coherence with the rest of the brain, equivalent to lower FCD) is replaced by the globally coherent PL state 1 (equivalent to higher FCD).

At the cellular level, studies have shown that 5-HT_2A receptors are highly expressed in layer 5 pyramidal neurons (Andrade and Weber, 2010). Stimulation of 5-HT_2A receptors induces a slow depolarization of pyramidal cells via G-protein coupled signaling pathways as well as the inhibition of the calcium-activated after-hyperpolarization (Andrade, 2011); together these changes make these neurons more likely to fire. Our results are therefore consistent with a scenario in which psilocybin increases the excitability of layer 5 pyramidal neurons broadly in the brain, wherever they are expressed. However, if the expression of 5-HT_2A receptors is especially high in lateral frontoparietal areas, then the dysregulated excitation that follows their stimulation (Carhart-Harris, 2019) may, in turn, destabilize this particular system, weakening its stability within the brain's broader repertoire. Supporting this idea, a recent modeling study showed that while the brain's structural connectome plays an important role in constraining and predicting brain dynamics, an additional layer of explanatory power is offered by considering energetic constraints on those dynamics (Gu et al., 2018). Because maintaining functional connections is energetically costly, the brain normally favors efficient neural codes and wiring patterns (Attwell and Laughlin, 2001; Bullmore and Sporns, 2012; Lord et al., 2013). In this framework, the metabolic costs of establishing and maintaining functional connections between anatomically distributed neurons in this Frontoparietal network may become too high following the psilocybin infusion, and brain activity might spend less time in this particular state as a result.

Consistent with this explanation, our analysis revealed increased stability of a functionally non-specific, globally coherent state (PL state 1) following the psilocybin injection. Previous work has indeed shown that the serotonergic psychedelics psilocybin and LSD induce an alternative type of functional integration, characterized by greater global integration (Petri et al., 2014; Roseman et al., 2014; Tagliazucchi et al., 2016) as well as increases in the metastability of brain states (Carhart-Harris et al., 2014). The results presented here not only align with these prior findings but also shed further light on their functional significance. Here we found that a globally coherent state, i.e., PL state 1, increased its predominance under psilocybin, not only statistically over time, but also dynamically, as revealed by the state-to-state switching matrix, where it became the 'state of choice' for inter-state transitions. In other words, under psilocybin the brain spent more time in a globally coherent, highly integrated state. This state therefore served as the dominant attractor, absorbing the majority of inter-state transitions. Residence in this particular, functionally non-specific state is consistent with a flattened attractor landscape where the dominant attractor reflects the relative absence of functionally specific attractors (Carhart-Harris and Friston, 2019). As was implied by previous analyses, this change in global brain dynamics may enable atypical patterns of interregional communication to occur in the psychedelic state (Petri et al., 2014) alongside the breakdown of within-network integrity in favor of globally integrated brain dynamics (Muthukumaraswamy et al., 2013; Carhart-Harris et al., 2014; Petri et al., 2014; Roseman et al., 2014; Tagliazucchi et al., 2016).

In addition to changing the relative stability of PL states and the directionality of state transitions, we found that psilocybin increased the global metastability of brain dynamics. This result is consistent with prior fMRI investigations of the same psilocybin dataset in which a greater diversity of functional connectivity motifs was observed after psilocybin infusion within a small network consisting of the bilateral hippocampal and anterior cingulate brain regions, thus reflecting increased variability in the collective repertoire of metastable states (Tagliazucchi et al., 2014). This prior analysis had purposely limited the number of regions included in order to perform an exhaustive counting of all possible functional connectivity states. Another study had found increased metastability in several canonical resting-state networks as determined by the variance in the networks' intrinsic synchrony over time (relationship between the mean and variance of the signal over time) (Carhart-Harris et al., 2014). Here we used a complementary analytical approach to demonstrate that increased metastability of BOLD phase coherence under psilocybin is generalizable at the whole-brain scale by calculating the standard deviation of the order parameter of the system over time.

The advantage of using a within-subject design (instead of a between-subject design) is that it improves the statistical power of the results. Indeed, despite the fact that there were only 9 subjects included in this fMRI dataset, our findings were strongly validated by the replication of results in the non-drug conditions across the two scanning days (i.e., no significant changes in the occupancy of network states in the placebo scanning session vs pre-psilocybin infusion), as well as the strong statistical effects observed in specific PL states after psilocybin injection (with p < 0.0005, surviving correction for the number of independent hypotheses tested). Importantly, since LEiDA compares the statistics between a discrete number of whole-brain PL patterns rather than individual functional links, the number of independent hypotheses tested is reduced from N(N-1)/2 (number of unique links in an NxN symmetric matrix) to the number of states K, which overcomes several of the statistical limitations associated with multiple comparisons in whole-brain dynamic connectivity analysis. Furthermore, given the instantaneous nature of LEiDA, we detected significant between-condition differences despite the fact that each scanning condition consisted only of 100 fMRI volumes, which would likely not have been sufficient for conventional sliding-window analysis (Hindriks et al., 2015). Moreover, we have demonstrated that the changes in the fractional occupancy of specific PL states after the psilocybin infusion was unrelated to participant motion during the scan.

The significant overlap between the BOLD phase-locking patterns detected automatically with LEiDA (from just 1800 fMRI volumes on 9 subjects) and the canonical resting-state networks obtained from ICA analysis on >1000 subjects (Yeo et al., 2011) suggests that the relative phase of the BOLD signals is sensitive to the same underlying patterns captured using ICA procedures, with the advantage that the networks can have spatial overlap and occur transiently at different frames in time. Although our approach did not reveal all the 7 RSNs reported in Yeo et al. (2011) - for instance, none of our PL states revealed the Dorsal Attention Network - this difference could be due to our limited sample size and the rather coarse parcellation scheme selected (AAL). Extending this methodological approach to parcellation schemes with a higher number of more homogenous areas (Craddock et al., 2013; Shen et al., 2013; Glasser et al., 2016) is likely to reveal more fine-grained PL states and deserves further exploration. Considering larger sample sizes and higher spatial resolution has the caveat of increasing the computational cost, but as long as the clustering algorithm adequately converges, it is likely to provide additional insights on the mechanisms underlying the formation of these BOLD phase-locking patterns and their relation to well-established RSNs. Importantly, using the cosine function, we focus on phase-locked oscillatory states and do not consider any directed causal influence between BOLD signals, which is fundamentally different to approaches considering asymmetric connectivity measures (Gilson et al., 2016). Regarding the temporal dynamics, we note that the LEiDA approach allows for the analysis of additional temporal properties of BOLD phase-locking patterns, including the dwell time of each state, as in Cabral et al. (2017b) and more recently in Figueroa et al. (2019). However, the temporal resolution used in the current dataset (TR = 3 s) was too slow to capture more rapid dynamics, which are suggested by MEG studies to operate at time scales of approximately 200 ms (Baker et al., 2014; Vidaurre et al., 2016). As such, we limited our study to the probabilities of occurrence and to the transition profiles occurring over these longer time scales.

The therapeutic potential of psilocybin in psychiatry has recently generated much interest. Indications of efficacy have been reported for conditions including treatment-resistant depression, anxiety related to end-of-life care and addictive disorders (Griffiths et al., 2006; Grob et al., 2011; Johnson et al., 2014; Carhart-Harris et al., 2016). The neural mechanisms underlying these clinical benefits, however, remain unclear. The present results provide some of the first evidence that a psychoactive compound modulates the brain's dynamical repertoire by selectively destabilizing a brain functional network (i.e., frontoparietal control system) and promoting transitions towards a globally coherent PL state. We believe this represents an exciting first step towards bridging molecular pharmacodynamics and dynamical measures of brain function in macro-scale networks, not only for serotonergic psychedelics but potentially across other classes of drugs and neuromodulators. This could have exciting implications for the design of novel therapeutics for neuropsychiatric disorders informed by patho-connectomics, which may be generally understood in terms of targeting the specific PL states affected by a particular disorder with neuromodulatory approaches. Moreover, despite the small sample size and the use of a coarse behavioral measure, we were nevertheless able to show strong associations between the fractional occupancy of certain PL states and the subjective intensity of the psychedelic experience. This interesting proof of concept study represents, to our knowledge, one of the first attempts at linking molecular pharmacodynamics, dynamical brain measures and behavior. We hope this work paves the way for fruitful future work in this direction.

In summary, the present study used a novel, data-driven dynamical functional connectivity analysis (LEiDA) to investigate how psilocybin rapidly modulates the exploration of the brain's repertoire of functional network states under task-free conditions. We found a PL state closely corresponding to a canonical frontoparietal control system to be markedly destabilized by the compound, while transitions toward a globally coherent PL state were enhanced. We also found an increase in the metastability of global brain dynamics following the psilocybin infusion. Taken together, these findings are consistent with prior neuroimaging studies suggesting that a different type of brain integration and increased neural signal complexity underlie the psychedelic state. The present results also suggest that characterizing neuromodulatory effects on the brain's dynamical repertoire may help guide circuit-specific pharmacological interventions in neuropsychiatric disorders.

Conflicts of interest
The authors declare no conflict of interest.

Acknowledgments

LDL is supported by the Canadian Institutes of Health Research, Canadian Centennial Scholarship Fund and the Mann Senior Scholarship from Hertford College, University of Oxford. JC is supported by the Portuguese Foundation for Science and Technology and by the project FronThera NORTE-01-0145-FEDER-000023 from the Northern Portugal Regional Operational Programme (NORTE 2020) under the Portugal 2020 Partnership Agreement through the European Regional Development Fund (FEDER). Dr. PE acknowledges support from the EPSRC award EP/N014529/1 funding the EPSRC Centre for Mathematics of Precision Healthcare at Imperial. MLK and SA are supported by the ERC Consolidator Grant CAREGIVING (n. 615539). MLK is also supported by the Center for Music in the Brain, funded by the Danish National Research Foundation (DNRF117). G.D. is supported by the Spanish Research Project PSI2016-75688-P (AEI/FEDER) and by the European Union's Horizon 2020 Framework Programme for Research and Innovation under the Specific Grant Agreement No. 785907 (Human Brain Project SGA2). RCH is supported by the Alex Mosley Charitable Trust, Ad Astra Chandaria Foundation and Tamas Foundation. The authors are grateful to the Beckley Foundation for funding the original research study, on which Amanda Feilding, Director of the Beckley Foundation, was a vital initiator and collaborative partner.

Appendix A. Supplementary data
Supplementary data to this article can be found online at https://doi.org/10.1016/j.neuroimage.2019.05.060.

References

1. Acebron, J.A., Bonilla, L.L., Vicente, C.J.P., Ritort, F., Spigler, R., 2005. The Kuramoto model: a simple paradigm for synchronization phenomena. Rev. Mod. Phys. 77, 137.
2. Acebron, J.A., Perales, A., Spigler, R., 2001. Bifurcations and global stability of synchronized stationary states in the Kuramoto model for oscillator populations. Phys. Rev. E, Stat. Nonlinear Soft Matter Phys. 64, 016218.
3. Andrade, R., 2011. Serotonergic regulation of neuronal excitability in the prefrontal cortex. Neuropharmacology 61, 382-386.
4. Andrade, R., Weber, E.T., 2010. Htr2a gene and 5-HT_2A receptor expression in the cerebral cortex studied using genetically modified mice. Front. Neurosci. 4, 36.
5. Atasoy, S., Vohryzek, J., Deco, G., Carhart-Harris, R.L., Kringelbach, M.L., 2018. Common neural signatures of psychedelics: frequency-specific energy changes and repertoire expansion revealed using connectome-harmonic decomposition. Prog. Brain Res. 242, 97-120.
6. Attwell, D., Laughlin, S.B., 2001. An energy budget for signaling in the grey matter of the brain. J. Cereb. Blood Flow Metab. 21, 1133-1145.
7. Baker, A.P., Brookes, M.J., Rezek, I.A., Smith, S.M., Behrens, T., Smith, P.J.P., Woolrich, M., 2014. Fast transient networks in spontaneous human brain activity. eLife 3.
8. Bayne, T., Carter, O., 2018. Dimensions of consciousness and the psychedelic state. Neurosci. Conscious. 2018 niy008.
9. Beckmann, C.F., DeLuca, M., Devlin, J.T., Smith, S.M., 2005. Investigations into resting-state connectivity using independent component analysis. Philos. Trans. R. Soc. Lond. Ser. B Biol. Sci. 360, 1001-1013.
10. Beckmann, C.F., Smith, S.M., 2004. Probabilistic independent component analysis for functional magnetic resonance imaging. IEEE Trans. Med. Imaging 23, 137-152.
11. Beliveau, V., Ganz, M., Feng, L., Ozenne, B., Højgaard, L., Fisher, P.M., Svarer, C., Greve, D.N., Knudsen, G.M., 2017. A high-resolution in vivo atlas of the human brain's serotonin system. J. Neurosci. 37, 120-128.
12. Biswal, B., Zerrin Yetkin, F., Haughton, V.M., Hyde, J.S., 1995. Functional connectivity in the motor cortex of resting human brain using echo-planar MRI. Magn. Reson. Med. 34, 537-541.
13. Bressler, S.L., Menon, V., 2010. Large-scale brain networks in cognition: emerging methods and principles. Trends Cognit. Sci. 14, 277-290.
14. Brookes, M.J., Hale, J.R., Zumer, J.M., Stevenson, C.M., Francis, S.T., Barnes, G.R., Owen, J.P., Morris, P.G., Nagarajan, S.S., 2011. Measuring functional connectivity using MEG: methodology and comparison with fcMRI. Neuroimage 56, 1082-1104.
15. Bullmore, E., Sporns, O., 2012. The economy of brain network organization. Nat. Rev. Neurosci. 13, 336.
16. Cabral, J., Hugues, E., Sporns, O., Deco, G., 2011. Role of local network oscillations in resting-state functional connectivity. Neuroimage 57, 130-139.
17. Cabral, J., Kringelbach, M., Deco, G., 2017a. Functional connectivity dynamically evolves on multiple time-scales over a static structural connectome: models and mechanisms. Neuroimage 160, 84-96.
18. Cabral, J., Kringelbach, M.L., Deco, G., 2014a. Exploring the network dynamics underlying brain activity during rest. Prog. Neurobiol. 114, 102-131.
19. Cabral, J., Luckhoo, H., Woolrich, M., Joensson, M., Mohseni, H., Baker, A., Kringelbach, M.L., Deco, G., 2014b. Exploring mechanisms of spontaneous functional connectivity in MEG: how delayed network interactions lead to structured amplitude envelopes of band-pass filtered oscillations. Neuroimage 90, 423-435.
20. Cabral, J., Vidaurre, D., Marques, P., Magalhaes, R., Silva Moreira, P., Miguel Soares, J., Deco, G., Sousa, N., Kringelbach, M.L., 2017b. Cognitive performance in healthy older adults relates to spontaneous switching between states of functional connectivity during rest. Sci. Rep. 7, 5135.
21. Carhart-Harris, R., Friston, K.J., 2019 July. REBUS and the anarchic brain: towards a unified model of the brain action of psychedelics. Pharmacol. Rev. 71, 1-28. https://doi.org/10.1124/pr.118.017160.
22. Carhart-Harris, R.L., 2019. How do psychedelics work? Curr. Opin. Psychiatr. 32, 16-21.
23. Carhart-Harris, R.L., Bolstridge, M., Rucker, J., Day, C.M., Erritzoe, D., Kaelen, M., Bloomfield, M., Rickard, J.A., Forbes, B., Feilding, A., 2016. Psilocybin with psychological support for treatment-resistant depression: an open-label feasibility study. The Lancet Psychiatr. 3, 619-627.
24. Carhart-Harris, R.L., Erritzoe, D., Williams, T., Stone, J.M., Reed, L.J., Colasanti, A., Tyacke, R.J., Leech, R., Malizia, A.L., Murphy, K., 2012. Neural correlates of the psychedelic state as determined by fMRI studies with psilocybin. Proc. Natl. Acad. Sci. Unit. States Am. 109, 2138-2143.
25. Carhart-Harris, R.L., Leech, R., Hellyer, P.J., Shanahan, M., Feilding, A., Tagliazucchi, E., Chialvo, D.R., Nutt, D., 2014. The entropic brain: a theory of conscious states informed by neuroimaging research with psychedelic drugs. Front. Hum. Neurosci. 8.
26. Cavanna, F., Vilas, M.G., Palmucci, M., Tagliazucchi, E., 2017. Dynamic functional connectivity and brain metastability during altered states of consciousness. Neuroimage 180, 383-395.
27. Christoff, K., Irving, Z.C., Fox, K.C., Spreng, R.N., Andrews-Hanna, J.R., 2016. Mind-wandering as spontaneous thought: a dynamic framework. Nat. Rev. Neurosci. 17, 718-731.
28. Cole, M.W., Schneider, W., 2007. The cognitive control network: integrated cortical regions with dissociable functions. Neuroimage 37, 343-360.
29. Craddock, R.C., Jbabdi, S., Yan, C.-G., Vogelstein, J., Castellanos, F.X., Di Martino, A., Kelly, C., Heberlein, K., Colcombe, S., Milham, M.P., 2013. Imaging human connectomes at the macroscale. Nat. Methods 10, 524-539.
30. Damoiseaux, J.S., Rombouts, S.A., Barkhof, F., Scheltens, P., Stam, C.J., Smith, S.M., Beckmann, C.F., 2006. Consistent resting-state networks across healthy subjects. Proc. Natl. Acad. Sci. U.S.A. 103, 13848-13853.
31. De Luca, M., Beckmann, C.F., De Stefano, N., Matthews, P.M., Smith, S.M., 2006. fMRI resting state networks define distinct modes of long-distance interactions in the human brain. Neuroimage 29, 1359-1367.
32. Deco, G., Jirsa, V., McIntosh, A.R., Sporns, O., Kotter, R., 2009. Key role of coupling, delay, and noise in resting brain fluctuations. Proc. Natl. Acad. Sci. U.S.A. 106, 10302-10307.
33. Deco, G., Kringelbach, M., 2016. Metastability and coherence: extending the communication through coherence hypothesis using a whole-brain computational perspective. Trends Neurosci. 39, 432.
34. Dosenbach, N.U., Fair, D.A., Cohen, A.L., Schlaggar, B.L., Petersen, S.E., 2008. A dual-networks architecture of top-down control. Trends Cognit. Sci. 12, 99-105.
35. Figueroa, C.A., Cabral, J., RJT, Mocking, Rapuano, K.M., van Hartevelt, T.J., Deco, G., Expert, P., Schene, A.H., Kringelbach, M.L., Ruhe, H.G., 2019. Altered ability to access a clinically relevant control network in patients remitted from major depressive disorder. Hum. Brain Mapp. 40, 2771-2786. https://doi.org/10.1002/hbm.24559.
36. Gilson, M., Moreno-Bote, R., Ponce-Alvarez, A., Ritter, P., Deco, G., 2016. Estimation of directed effective connectivity from fMRI functional connectivity hints at asymmetries of cortical connectome. PLoS Comput. Biol. 12, e1004762.
37. Glasser, M.F., Coalson, T.S., Robinson, E.C., Hacker, C.D., Harwell, J., Yacoub, E., Ugurbil, K., Andersson, J., Beckmann, C.F., Jenkinson, M., Smith, S.M., Van Essen, D.C., 2016. A multi-modal parcellation of human cerebral cortex. Nature 536, 171-178.
38. Glerean, E., Salmi, J., Lahnakoski, J.M., Jaaskelainen, I.P., Sams, M., 2012. Functional magnetic resonance imaging phase synchronization as a measure of dynamic functional connectivity. Brain Connect. 2, 91-101.
39. González-Maeso, J., Weisstaub, N.V., Zhou, M., Chan, P., Ivic, L., Ang, R., Lira, A., Bradley-Moore, M., Ge, Y., Zhou, Q., 2007. Hallucinogens recruit specific cortical 5-HT 2A receptor-mediated signaling pathways to affect behavior. Neuron 53, 439-452.
40. Griffiths, R.R., Richards, W.A., McCann, U., Jesse, R., 2006. Psilocybin can occasion mystical-type experiences having substantial and sustained personal meaning and spiritual significance. Psychopharmacology 187, 268-283.
41. Grob, C.S., Danforth, A.L., Chopra, G.S., Hagerty, M., McKay, C.R., Halberstadt, A.L., Greer, G.R., 2011. Pilot study of psilocybin treatment for anxiety in patients with advanced-stage cancer. Arch. Gen. Psychiatr. 68, 71-78.
42. Gu, S., Cieslak, M., Baird, B., Muldoon, S.F., Grafton, S.T., Pasqualetti, F., Bassett, D.S., 2018. The energy landscape of neurophysiological activity implicit in brain network structure. Sci. Rep. 8, 2507.
43. Gutierrez-Barragan, D., Basson, M.A., Panzeri, S., Gozzi, A., 2018. Oscillatory Brain States Govern Spontaneous fMRI Network Dynamics bioRxiv 393389.
44. Hansen, E.C., Battaglia, D., Spiegler, A., Deco, G., Jirsa, V.K., 2015. Functional connectivity dynamics: modeling the switching behavior of the resting state. Neuroimage 105, 525-535.
45. Hellyer, P.J., Scott, G., Shanahan, M., Sharp, D.J., Leech, R., 2015. Cognitive flexibility through metastable neural dynamics is disrupted by damage to the structural connectome. J. Neurosci.: Off. J. Soc. Neurosci. 35, 9050-9063.
46. Hindriks, R., Adhikari, M.H., Murayama, Y., Ganzetti, M., Mantini, D., Logothetis, N.K., Deco, G., 2015. Can sliding-window correlations reveal dynamic functional connectivity in resting-state fMRI? Neuroimage 127, 242-256.
47. Hipp, J.F., Hawellek, D.J., Corbetta, M., Siegel, M., Engel, A.K., 2012. Large-scale cortical correlation structure of spontaneous oscillatory activity. Nat. Neurosci. 15, 884-890.
48. Jenkinson, M., Bannister, P., Brady, M., Smith, S., 2002. Improved optimization for the robust and accurate linear registration and motion correction of brain images. Neuroimage 17, 825-841.
49. Johnson, M.W., Garcia-Romeu, A., Cosimano, M.P., Griffiths, R.R., 2014. Pilot study of the 5-HT_2AR agonist psilocybin in the treatment of tobacco addiction. J. Psychopharmacol. 28, 983-992.
50. Karahanoǧlu, F.I., Van De Ville, D., 2015. Transient brain activity disentangles fMRI resting-state dynamics in terms of spatially and temporally overlapping networks. Nat. Commun. 6, 7751.
51. Kometer, M., Schmidt, A., Jäncke, L., Vollenweider, F.X., 2013. Activation of serotonin 2A receptors underlies the psilocybin-induced effects on α oscillations, N170 visual-evoked potentials, and visual hallucinations. J. Neurosci. 33, 10544-10551.
52. Lempel, A., Ziv, J., 1976. On the complexity of finite sequences. IEEE Trans. Inf. Theory 22, 75-81.
53. Liu, X., Duyn, J.H., 2013. Time-varying functional network information extracted from brief instances of spontaneous brain activity. Proc. Natl. Acad. Sci. U.S.A. 110, 4392-4397.
54. Lord, L.-D., Expert, P., Huckins, J.F., Turkheimer, F.E., 2013. Cerebral energy metabolism and the brain's functional network architecture: an integrative review. J. Cereb. Blood Flow Metab. 33, 1347-1354.
55. Lord, L.-D., Stevner, A.B., Deco, G., Kringelbach, M.L., 2017. Understanding principles of integration and segregation using whole-brain computational connectomics: implications for neuropsychiatric disorders. Phil. Trans. R Soc. A 375, 20160283.
56. Mantini, D., Perrucci, M.G., Del Gratta, C., Romani, G.L., Corbetta, M., 2007. Electrophysiological signatures of resting state networks in the human brain. Proc. Natl. Acad. Sci. U.S.A. 104, 13170-13175.
57. McKenna, D.J., Repke, D.B., Lo, L., Peroutka, S.J., 1990. Differential interactions of indolealkylamines with 5-hydroxytryptamine receptor subtypes. Neuropharmacology 29, 193-198.
58. Meltzer, H., Massey, B., Horiguchi, M., 2012. Serotonin receptors as targets for drugs useful to treat psychosis and cognitive impairment in schizophrenia. Curr. Pharmaceut. Biotechnol. 13, 1572-1586.
59. Menon, V., 2011. Large-scale brain networks and psychopathology: a unifying triple network model. Trends Cognit. Sci. 15, 483-506.
60. Musso, F., Brinkmeyer, J., Mobascher, A., Warbrick, T., Winterer, G., 2010. Spontaneous brain activity and EEG microstates. A novel EEG/fMRI analysis approach to explore resting-state networks. Neuroimage 52, 1149-1161.
61. Muthukumaraswamy, S.D., Carhart-Harris, R.L., Moran, R.J., Brookes, M.J., Williams, T.M., Errtizoe, D., Sessa, B., Papadopoulos, A., Bolstridge, M., Singh, K.D., 2013. Broadband cortical desynchronization underlies the human psychedelic state. J. Neurosci. 33, 15171-15183.
62. Newman, M.E., 2006. Finding community structure in networks using the eigenvectors of matrices. Phys. Rev. E, Stat. Nonlinear Soft Matter Phys. 74, 036104.
63. Passie, T., Seifert, J., Schneider, U., Emrich, H.M., 2002. The pharmacology of psilocybin. Addict. Biol. 7, 357-364.
64. Petri, G., Expert, P., Turkheimer, F., Carhart-Harris, R., Nutt, D., Hellyer, P.J., Vaccarino, F., 2014. Homological scaffolds of brain functional networks. J. R. Soc. Interface 11, 20140873.
65. Ponce-Alvarez, A., Deco, G., Hagmann, P., Romani, G.L., Mantini, D., Corbetta, M., 2015. Resting-state temporal synchronization networks emerge from connectivity topology and heterogeneity. PLoS Comput. Biol. 11, e1004100.
66. Roberts, J.A., Gollo, L.L., Abeysuriya, R.G., Roberts, G., Mitchell, P.B., Woolrich, M.W., Breakspear, M., 2019. Metastable brain waves. Nat. Commun. 10, 1056.
67. Roseman, L., Leech, R., Feilding, A., Nutt, D.J., Carhart-Harris, R.L., 2014. The effects of psilocybin and MDMA on between-network resting state functional connectivity in healthy volunteers. Front. Hum. Neurosci. 8, 204.
68. Schartner, M.M., Carhart-Harris, R.L., Barrett, A.B., Seth, A.K., Muthukumaraswamy, S.D., 2017. Increased spontaneous MEG signal diversity for psychoactive doses of ketamine, LSD and psilocybin. Sci. Rep. 7, 46421.
69. Seeley, W.W., Menon, V., Schatzberg, A.F., Keller, J., Glover, G.H., Kenna, H., Reiss, A.L., Greicius, M.D., 2007. Dissociable intrinsic connectivity networks for salience processing and executive control. J. Neurosci.: Off. J. Soc. Neurosci. 27, 2349-2356.
70. Shanahan, M., 2010. Metastable chimera states in community-structured oscillator networks. Chaos 20, 013108.
71. Shen, X., Tokoglu, F., Papademetris, X., Constable, R.T., 2013. Groupwise whole-brain parcellation from resting-state fMRI data for network node identification. Neuroimage 82, 403-415.
72. Smallwood, J., Brown, K., Baird, B., Schooler, J.W., 2012. Cooperation between the default mode network and the frontal-parietal network in the production of an internal train of thought. Brain Res. 1428, 60-70.
73. Smith, S.M., 2002. Fast robust automated brain extraction. Hum. Brain Mapp. 17, 143-155.
74. Tagliazucchi, E., Balenzuela, P., Fraiman, D., Chialvo, D.R., 2012. Criticality in large-scale brain FMRI dynamics unveiled by a novel point process analysis. Front. Physiol. 3, 15.
75. Tagliazucchi, E., Carhart-Harris, R., Leech, R., Nutt, D., Chialvo, D.R., 2014. Enhanced repertoire of brain dynamical states during the psychedelic experience. Hum. Brain Mapp. 35, 5442-5456.
76. Tagliazucchi, E., Roseman, L., Kaelen, M., Orban, C., Muthukumaraswamy, S.D., Murphy, K., Laufs, H., Leech, R., McGonigle, J., Crossley, N., 2016. Increased global functional connectivity correlates with LSD-induced ego dissolution. Curr. Biol. 26, 1043-1050.
77. Tewarie, P., Abeysuriya, R., Byrne, Á., O'Neill, G.C., Sotiropoulos, S.N., Brookes, M.J., Coombes, S., 2019. How do spatially distinct frequency specific MEG networks emerge from one underlying structural connectome? The role of the structural eigenmodes. Neuroimage 186, 211-220.
78. Tononi, G., Edelman, G.M., 1998. Consciousness and complexity. Science 282, 1846-1851.
79. Turton, S., Nutt, D., Carhart-Harris, R., 2014. A qualitative report on the subjective experience of intravenous psilocybin administered in an FMRI environment. Curr. Drug Abuse Rev. 7, 117-127.
80. Tzourio-Mazoyer, N., Landeau, B., Papathanassiou, D., Crivello, F., Etard, O., Delcroix, N., Mazoyer, B., Joliot, M., 2002. Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain. Neuroimage 15, 273-289.
81. Varela, F.J., 1979. Principles of Biological Autonomy.
82. Vidaurre, D., Quinn, A.J., Baker, A.P., Dupret, D., Tejero-Cantero, A., Woolrich, M.W., 2016. Spectrally resolved fast transient brain states in electrophysiological data. Neuroimage 126, 81-95.
83. Vidaurre, D., Smith, S.M., Woolrich, M.W., 2017. Brain network dynamics are hierarchically organized in time. Proc. Natl. Acad. Sci. Unit. States Am. 114, 12827-12832.
84. Vincent, J.L., Kahn, I., Snyder, A.Z., Raichle, M.E., Buckner, R.L., 2008. Evidence for a frontoparietal control system revealed by intrinsic functional connectivity. J. Neurophysiol. 100, 3328-3342.
85. Vollenweider, F.X., Vollenweider-Scherpenhuyzen, M.F., Bäbler, A., Vogel, H., Hell, D., 1998. Psilocybin induces schizophrenia-like psychosis in humans via a serotonin-2 agonist action. Neuroreport 9, 3897-3902.
86. Watanabe, T., Hirose, S., Wada, H., Imai, Y., Machida, T., Shirouzu, I., Konishi, S., Miyashita, Y., Masuda, N., 2014. Energy landscapes of resting-state brain networks. Front. Neuroinf. 8, 12.
87. Xin, F., Lei, X., 2015. Competition between frontoparietal control and default networks supports social working memory and empathy. Soc. Cognit. Affect Neurosci. 10, 1144-1152.
88. Yeo, B.T., Krienen, F.M., Sepulcre, J., Sabuncu, M.R., Lashkari, D., Hollinshead, M., Roffman, J.L., Smoller, J.W., Zöllei, L., Polimeni, J.R., 2011. The organization of the human cerebral cortex estimated by intrinsic functional connectivity. J. Neurophysiol. 106, 1125.