# Dynamic coupling of whole-brain neuronal and neurotransmitter systems

^{a}Department of Psychiatry, University of Oxford, Oxford OX3 7JX, United Kingdom;^{b}Center for Music in the Brain, Department of Clinical Medicine, Aarhus University, DK-8000 Aarhus C, Denmark;^{c}Life and Health Sciences Research Institute, School of Medicine, University of Minho, 4710-057 Braga, Portugal;^{d}Centre for Eudaimonia and Human Flourishing, University of Oxford, Oxford OX1 2JD, United Kingdom;^{e}Center for Brain and Cognition, Computational Neuroscience Group, Universitat Pompeu Fabra, 08018 Barcelona, Spain;^{f}Department of Information and Communication Technologies, Universitat Pompeu Fabra, 08018 Barcelona, Spain;^{g}Neurobiology Research Unit, Rigshospitalet, DK-2100 Copenhagen, Denmark;^{h}Center for Integrated Molecular Brain Imaging, Rigshospitalet, DK-2100 Copenhagen, Denmark;^{i}Faculty of Health and Medical Sciences, Copenhagen University, DK-2100 Copenhagen, Denmark;^{j}The Centre for Psychedelic Research, Department of Brain Sciences, Imperial College London, London W12 0NN, United Kingdom;^{k}Semel Institute for Neuroscience and Human Behavior, University of California, Los Angeles, CA 90024;^{l}Max Planck Institute for Biological Cybernetics, 72076 Tübingen, Germany;^{m}Imaging Science and Biomedical Engineering, University of Manchester, Manchester M13 9PT, United Kingdom;^{n}Institució Catalana de la Recerca i Estudis Avançats, 08010 Barcelona, Spain;^{o}Department of Neuropsychology, Max Planck Institute for Human Cognitive and Brain Sciences, 04103 Leipzig, Germany;^{p}School of Psychological Sciences, Monash University, Clayton, Melbourne, VIC 3800, Australia

See allHide authors and affiliations

Contributed by Nikos K. Logothetis, February 23, 2020 (sent for review December 10, 2019; reviewed by Rodrigo Quian Quiroga and Olaf Sporns)

## Significance

In a technical tour de force, we have created a framework demonstrating the underlying fundamental principles of bidirectional coupling of neuronal and neurotransmitter dynamical systems. Specifically, in the present study, we combined multimodal neuroimaging data to causally explain the functional effects of specific serotoninergic receptor (5-HT_{2A}R) stimulation with psilocybin in healthy humans. Longer term, this could provide a better understanding of why psilocybin is showing considerable promise as a therapeutic intervention for neuropsychiatric disorders including depression, anxiety, and addiction.

## Abstract

Remarkable progress has come from whole-brain models linking anatomy and function. Paradoxically, it is not clear how a neuronal dynamical system running in the fixed human anatomical connectome can give rise to the rich changes in the functional repertoire associated with human brain function, which is impossible to explain through long-term plasticity. Neuromodulation evolved to allow for such flexibility by dynamically updating the effectivity of the fixed anatomical connectivity. Here, we introduce a theoretical framework modeling the dynamical mutual coupling between the neuronal and neurotransmitter systems. We demonstrate that this framework is crucial to advance our understanding of whole-brain dynamics by bidirectional coupling of the two systems through combining multimodal neuroimaging data (diffusion magnetic resonance imaging [dMRI], functional magnetic resonance imaging [fMRI], and positron electron tomography [PET]) to explain the functional effects of specific serotoninergic receptor (5-HT_{2A}R) stimulation with psilocybin in healthy humans. This advance provides an understanding of why psilocybin is showing considerable promise as a therapeutic intervention for neuropsychiatric disorders including depression, anxiety, and addiction. Overall, these insights demonstrate that the whole-brain mutual coupling between the neuronal and the neurotransmission systems is essential for understanding the remarkable flexibility of human brain function despite having to rely on fixed anatomical connectivity.

Human connectomics has been very successful in revealing how function arises from structure (1, 2), and by showing how anatomy can give rise to a complex dynamical neuronal system as measured with multimodal neuroimaging (3⇓–5). Despite the attractiveness of the idea that function is shaped by anatomy, it is also clear that the single fixed structure of the anatomical connectome should not be able to give rise to the full palette and complexity of brain function. However, evolution has found a solution to this apparent paradox by dynamically modulating the connectome over time through neuromodulatory systems, enabling the richness of behaviors needed for survival. Indeed, the necessary dynamics of human brain function can be obtained by modulating the effective connectivity of the coupling over time, as proposed by Friston and many others (6, 7). Still, a principled and mechanistic description of the dynamic connectome must bring together the anatomical, neuronal, and neurotransmitter systems at the whole-brain level (8).

Here, we show how the mutual coupling between the neuronal and neurotransmitter systems is fundamental to understanding the dynamic connectome. This can be achieved through whole-brain modeling of multimodal neuroimaging data using a mutually coupled neuronal–neurotransmitter whole-brain model where the structural anatomical connectivity can be measured using diffusion magnetic resonance imaging (dMRI), the functional connectivity with functional magnetic resonance imaging (fMRI), and neurotransmission (receptor density) with positron electron tomography (PET). In this model, the synaptic/neuronal activity and the neurotransmitter diffusive system are portrayed in a realistic biophysical way by a set of separate dynamical equations, which are mutually coupled through the receptor maps and synaptic dynamics (neurotransmitters to neuronal), and the excitation of the projections from the brain regions producing the neurotransmitters (neuronal to neurotransmitters). The explicit linkage between this dual-coupled dynamical system yields a deeper understanding of the crucial mutual coupling between neuronal and neurotransmitters systems at the whole-brain level. We perceive this as a significant improvement compared to our previous unique whole-brain study, which had only one neuronal dynamical system influenced by static neurotransmitter concentrations modulating the neuronal gain (9). Specifically, we demonstrate the explanatory and predictive power of this mutually coupled whole-brain model by investigating the effects of psychedelics on brain activity.

Psilocybin, the prodrug of psilocin (4-OH-dimethyltryptamine), is a good model system to demonstrate the power of a mutually coupled whole-brain model, since it has been shown to act mainly through the serotonin 2A receptor (5-HT_{2A}R) (10), rather than more complex interactions between many receptors. The serotonin system works through the projections of the raphe nucleus. In addition, conveniently, the 5-HT receptor density maps have recently been mapped with PET (11). Here, we were interested in revealing the effects of mutual coupling of both neuronal and neurotransmitter systems on brain repertoire and specifically the effects of psilocybin on resting-state activity on healthy human participants.

Specifically, the bidirectional coupling of the neuronal and neurotransmitter systems was modeled in the following way: For the placebo condition, we used a standard whole-brain model to simulate the neuronal system, i.e., modeling spontaneous brain activity at the whole-brain level (measured with blood oxygen level-dependent [BOLD] fMRI), where each node represents a brain area and the links between them are represented by white matter connections (measured with dMRI). For the psilocybin condition, we mutually coupled the whole-brain neuronal and neurotransmitter systems by including an explicit description of the neurotransmitter dynamical system and the mutual coupling with the neuronal system. This was done by modeling the dynamics of the neurotransmitter system through simulating the release-and-reuptake dynamics, where the serotonin receptor density of each brain area is measured with PET. The neurotransmitter dynamics are then in turn coupled with the neuronal activity through the firing rate activity of the raphe nucleus, source of the serotonin neurotransmitter.

The relevance of the whole-brain modeling developed here is strongly supported by recent studies that have started to demonstrate the functional neuroanatomy underlying the experience of unconstrained cognition and enhanced mind-wandering reported following psilocybin (12⇓⇓–15). Due to its therapeutic action for the treatment of neuropsychiatric disorders such as depression, anxiety, and addiction (16⇓–18), psilocybin has begun to elicit significant interest from the neuropsychiatric community as a potential treatment (19). Long term, the *in silico* framework presented here for investigating the dynamical connectome has the potential to bring insights and help design interventions for brain disease including neuropsychiatric disorders, which are otherwise difficult to study with traditional animal models.

## Results

The main aim of the mutually coupled neuronal–neurotransmitter whole-brain model (shown in Fig. 1) is to investigate the tight interactions between these two different but mutually coupled dynamical whole-brain systems (Fig. 1*A*), which are fitted (Fig. 1*B*) to the empirical neuroimaging BOLD data (Fig. 1*C*). In other words, neuronal-only models that fit neuroimaging data using local node dynamics constrained by anatomical structural connectivity (20⇓⇓–23) are now described by a balanced dynamic mean field model that expresses consistently the time evolution of the ensemble activity of the different neural populations building up the spiking network (24, 25). The spontaneous activity of each single brain region is described by a network consisting of two pools of excitatory and inhibitory neurons (Fig. 1*D* and Eqs. **2**–**7** in *Materials and Methods*). The neurotransmitter system, on the other hand, uses a separate set of differential equations describing the dynamics of the neurotransmitter concentration level, given by the well-known Michaelis–Menten release-and-reuptake dynamics (Eq. **13** in *Materials and Methods* and shown in Fig. 1*E*) (26⇓–28). We couple the neurotransmitter and neuronal dynamical systems by means of the anatomical connectivity between the raphe nucleus and the rest of the brain, estimated using Human Connectome Project (HCP) dMRI (*Materials and Methods*). The explicit coupling between the neurotransmitter and the neuronal system is shown in Fig. 1*F* (and in Eqs. **14** and **15** in *Materials and Methods*). As can be seen, the neurotransmitter currents are applied to each region’s excitatory and inhibitory pools of neurons using the effectivity/conductivity parameters (*W*_{E} and *W*_{I}, respectively). In each region, the neurotransmitter currents are also scaled by each region’s receptor density (measured in vivo using PET). The reverse coupling from the neuronal to the neurotransmitter system is given by inserting in the Michaelis–Menten release-and-reuptake equation the neuronal population firing rate of the source brain region of the neurotransmitter spread from the raphe nucleus.

As a demonstration of the power of this general framework, we show how it can used to explain the specific case of modulation of the serotonin system by psilocybin by analyzing the neuroimaging data from a group of healthy participants being administered psilocybin i.v. (see details in *Materials and Methods*). This allowed us to study the relevant role of the coupling between the neuronal and neurotransmitter systems to understand the specificity of the 5-HT_{2A} receptor densities [mapped with PET in vivo (11)]. Specifically, the main source of serotonin neurotransmission is the raphe nucleus, which is essential for the correct coupling between both dynamical systems. We therefore used the structural connectivity from this region estimated using dMRI tractography (*Materials and Methods*).

We fitted the mutually coupled whole-brain model using our framework of describing a brain state as an ensemble or probabilistic “cloud” in a given state space (29). This cloud is of course not clustered into discrete states (30), but it has been shown that clustering can nevertheless be useful in providing so-called “metastable substates” that can significantly distinguish between brain states (31, 32). A brain state is determined by a collection of metastable substates, i.e., of time-varying pseudo-states resulting from the clustering process (33⇓–35).

We extracted this probabilistic metastable substate (PMS) space for the empirical psilocybin data (placebo and active condition) using the leading eigenvector dynamics analysis (LEiDA) method (described in detail in *Materials and Methods* and schematized in Fig. 2) (31). Briefly, we define a time-resolved dynamic FC (dFC) matrix by using BOLD phase coherence connectivity. In order to reduce the dimensionality of the problem, we compute the corresponding time-resolved leading eigenvector, which captures the dominant connectivity pattern of dFC(*t*) at time *t*. Then we detect a discrete number of metastable substates by clustering the dFC(*t*) across time points and subjects. The obtained *k* cluster centroids define the PMS space, for which we compute the probabilities, lifetimes, and transition probability between them for both the placebo and active conditions of psilocybin. The placebo and the active conditions of the psilocybin can be significantly distinguished by three substates. Fig. 3*A* shows three different substates of the PMS methodology with associated probabilities and transition probabilities between them for the placebo and active conditions of the psilocybin. As can be seen in the subplots, two substates (1 and 3) are significantly different between the two conditions (*P* < 10^{−4}) in terms of probability and substate 3 is significantly different for lifetimes (*P* < 10^{−2}). This demonstrates that the clustering approach is indeed useful for distinguishing brain states.

We first fitted the whole-brain model to the PMS space of the placebo condition of psilocybin using only the neuronal system (and thus without coupling the neurotransmitter system). We did this by fitting the minimum of the symmetrized Kullback–Leibler distance (KLD) between the empirical placebo condition PMS and modeled PMS (*Materials and Methods*). This yielded a global neuronal coupling parameter *G* = 1.6. Equally, we measured the Markov entropy distance (ME) (*Materials and Methods*). The optimal values were found to be KLD = 0.002 and ME = 0.05.

Second, we proceeded to study the role of the neurotransmitter dynamical system by coupling both systems in the model. This allowed us to predict the changes in brain dynamics under the effects of psilocybin. Specifically, the mutually coupled whole-brain model used the same centroids found in the neuroimaging empirical data and a two-dimensional (2D) matrix of coupling parameters *B*). We searched this matrix for the optimal fitting (global minimum) to the active condition of the psilocybin data by using the optimal symmetrized KLD between the empirical active psilocybin condition PMS and the modeled PMS. Nevertheless, it is clear from the deep diagonal ridge that there is an invariant ratio between *B* shows the matrices of the KLD and the error lifetimes of the substates in the empirical active psilocybin condition PMS and the modeled PMS as a function of coupling parameters *Materials and Methods*). The optimal description of the active psilocybin condition is found at the minimum of the KLD at *C* shows the PMS spaces for the uncoupled system (

At this optimal fit, the mutually coupled whole-brain model allows us to obtain further insights into the underlying dynamics of neurotransmission involved in psilocybin (in this case for serotonin). Fig. 4*A* shows a significant difference between the optimal fit and the uncoupled system in KLD (*P* < 10^{−6}). This clearly demonstrates the significance of coupling the neuronal and neurotransmitter systems. Fig. 4*B* further dissects this finding by showing a significant difference between the optimal fit and the optimal fit but where we have frozen the feedback dynamics from the neurotransmitter system to the neuronal system (*P* < 10^{−6}). This was done by allowing the coupling until the steady state is achieved, and then we take the average of the neurotransmitter variables and freeze these and cancel the feedback dynamics. With this uncoupling, we consider a simpler coupled system, where the neuronal dynamics are preserved but unable to influence the neurotransmitter dynamics that are frozen.

Fig. 4*C* shows the significant difference between using the empirical 5-HT_{2A} receptor densities across the regions at the optimal fit compared with the results of randomly shuffling the receptor densities (*P* < 10^{−4}). Similarly, as shown in Fig. 5*A*, we also tested the specificity of the receptor binding maps by comparing (at the optimal fit) the 5-HT_{2A} receptor with the other serotonin receptors, namely 5-HT_{1A}, 5-HT_{1B}, and 5-HT_{4}, which showed significant differences between 5-HT_{2A} and 5-HT_{1A} (*P* < 10^{−6}), 5-HT_{1B} (*P* < 10^{−6}), and 5-HT_{4} (*P* < 10^{−6}). This strongly confirms the main role this receptor plays in the effects of psilocybin (10). Finally, in Fig. 5*B*, we show the comparison to the 5-HTT, which is also significant (*P* < 0.02). 5-HTT is not a receptor but has been shown to play an important role for the treatment of depression.

## Discussion

Here, we have shown how a dynamic mutually coupled whole-brain model can address the major challenge in neuroscience, which is to explain the paradoxical flexibility of human brain function despite the reliance on a fixed anatomical connectome. One of the most important elements in this flexibility is the way that the fixed connectome can be modulated by neuromodulators to selectively change the balance of the excitation and inhibition of individual brain regions. Here, we modeled the dynamical mutual coupling between neuronal and neuromodulator systems at the whole-brain level. In particular, we implemented a mutually coupled dynamical system, which couples at the whole-brain level the neuronal system, which was modeled using a balanced dynamic mean field model (24, 25) and the neuromodulator system describing the dynamics of the neurotransmitter concentration levels (measured in vivo using PET) and modeled by the well-known Michaelis–Menten release-and-reuptake dynamics (26⇓–28). Here, as proof of principle, we consider the effects of psilocybin on the serotonin system and therefore use anatomical connectivity between the raphe nucleus and the rest of the brain to couple the two systems.

Overall, the results show that the interaction between these dynamical systems is fundamental for explaining the empirical data. In other words, the dynamic mutual interaction between neuronal and neuromodulator systems at the whole-brain level is important to fully explain the functional modulation of brain activity by psilocybin, a powerful psychedelic drug, acting on the serotonin system. This is especially important given the demonstrated ability of psilocybin to rebalance the human brain in treatment-resistant depression. The results provide evidence for how the integration of dMRI (anatomy), fMRI (functional neuronal activity), and PET (neurotransmitter system) at the whole-brain level is necessary for predicting properly brain dynamics as a result of the mutual coupling between a dual dynamical system. This expands on the rich existing experimental and theoretical research on the local effects of neurotransmitters (e.g., refs. 36⇓⇓⇓–40).

Specifically, the results provide insights into the underlying dynamics of neurotransmission involved in psilocybin. In terms of dynamics, we first uncoupled the neuromodulators from the neuronal systems and found that this produced a highly significant breakdown in fitting the empirical data (Fig. 4*A*). We then ran further simulations to investigate the role of specific parts of the dynamic coupling. The full mutually coupled whole-brain model ran until a steady state was achieved, upon which we then froze the feedback dynamics from neuromodulators to neuronal system by removing coupling through the raphe nucleus. This again produced a highly significant breakdown in fitting the empirical data (Fig. 4*B*).

In further sets of experiments designed to investigate the importance of the receptor distribution, we changed the distribution of regional receptor densities by 1) randomly shuffling the 5-HT_{2A} (Fig. 4*C*) and 2) replacing them with those of other serotonin receptors, known to be much less sensitive to psilocybin, namely 5-HT_{1A}, 5-HT_{1B}, and 5-HT_{4}. Again, this produced a highly significant breakdown in the ability to explain the empirical data (Fig. 5*B*). This result confirms the crucial, causative role of the precise anatomical location and density of 5-HT_{2A}.

Interestingly, as mentioned above, psilocybin has been demonstrated to play a role in rebalancing the human brain in treatment-resistant depression (19). The therapeutic action of psilocybin in depression is thought to depend on activation of serotonin 2A receptors, thereby initiating a multilevel plasticity (41). This is different from selective serotonin reuptake inhibitors, which are the most frequently prescribed antidepressant drug class and whose therapeutic action is thought to begin with reuptake blockade at the 5-HTT; psilocybin has no appreciable affinity or action at the 5-HTT. We were interested in further investigating this by replacing the receptor densities with the 5-HTT densities. We found that the fit with 5-HTT was significantly worse than 5-HT_{2A} (Fig. 5*C*) and supports the potential for psilocybin in rebalancing brain function in treatment-resistant depression (19).

More broadly, the mutually coupled whole-brain model and results shed important light on our understanding of human brain function. Over the last decade, the discovery of whole-brain modeling has led to important new insights for explaining the principles of human brain function based on human neuroimaging data. The initial breakthrough was the recognition that the anatomical connectivity of the human brain shapes function at the whole-brain level (1, 2). However, it also became clear relatively quickly that this was a by-product of the close link between anatomy and function, given that it is possible to describe static spontaneous activity (i.e., grand average functional spatial connectivity over longer periods) even in the absence of neuronal dynamics, i.e., combining the underlying anatomy with noise (42). Still, it was shown by further investigations that capturing the dynamics of the spatiotemporal whole-brain activity requires more sophisticated dynamical modeling (3, 35, 43). All of these whole-brain studies were initially solely based on neuronal modeling, but recently it was shown that even better results can be obtained when the dynamics of the neuronal system (through the neuronal gain of excitatory neurons) is influenced by static neuromodulation (9).

Here, we have demonstrated the importance of having a biophysically realistic dynamic neuromodulator system, and more importantly the need for full mutual coupling between the full dynamics of both neuronal and neuromodulation systems. In other words, this framework for mutually coupled whole-brain modeling involves the underlying anatomical connectivity and two mutually interacting dynamical neuronal and neuromodulator systems.

Overall, the framework put forward here is likely to be essential for accurately modeling and explaining mechanisms of human brain function in health and disease. The framework combines multimodal neuroimaging from dMRI (anatomy), fMRI (functional neuronal activity), and PET (neurotransmitter system) at the whole-brain level. This offers unique opportunities for investigating how to rebalance human brain activity in disease through pharmacological manipulation that can be discovered and tested *in silico* in the mutually coupled whole-brain model proposed here. This whole-brain modeling can thus be used to make predictions that can be tested causally with electromagnetic stimulation (deep brain stimulation [DBS] and transcranial magnetic stimulation) or pharmacological interventions. In the future, this framework could further extended by measuring even faster whole-brain dynamics (measured with magnetoencephalography) (44) and incorporating gene expression at the whole-brain level, offering even more exciting avenues for discovering new and effective ways of rebalancing the human brain in disease.

## Materials and Methods

We provide here the details on the pipeline used to integrate structural and functional connectivity (dMRI and fMRI) with neurotransmission (PET) in a model of the placebo and psilocybin response in healthy participants (summarized in Fig. 1). Please note that all structural, functional and neuromodulation data are integrated into the Automated Anatomical Labeling (AAL) parcellation:

1) Structural connectivity: probabilistic tractography derived from the dMRI;

2) Functional dynamics: functional dynamics described as probability metastable substates were estimated from the fMRI data obtained in the placebo and psilocybin condition;

3) Neurotransmitter receptor density: estimation of the density of the 5HT-2A receptors that has been obtained using PET;

4) Whole-brain neuronal model: description of the (uncoupled) neuronal dynamic mean field model used to fit the placebo condition;

5) Mutually coupled neuronal and neurotransmission whole-brain model: dynamic coupling of the neuronal and neurotransmission models to dynamically fit psilocybin condition;

6) Empirical fitting of mutually coupled whole-brain model to neuroimaging data: measuring the fit of model to the empirical neuroimaging data.

### Parcellation.

Based on our previous whole-brain studies, we used the AAL atlas but considered only the 90 cortical and subcortical noncerebellar brain regions (45). All structural, functional, and neuromodulation data were integrated using this atlas.

#### Structural connectivity.

##### Ethics.

The study was approved by the Research Ethics Committee of the Central Denmark Region (De Videnskabsetiske Komitéer for Region Midtjylland). Written informed consent was obtained from all participants prior to participation.

##### Participants and acquisition.

The dMRI data used in the study was collected from 16 healthy right-handed participants at Aarhus University, Denmark (11 men; mean age, 24.8 ± 2.5) (46, 47). We recorded the imaging data in a single session on a 3-T Siemens Skyra scanner at Center of Functionally Integrative Neuroscience, Aarhus University, Denmark. The structural MRI T1 scan used the following parameters: reconstructed matrix size, 256 × 256; voxel size of 1 mm^{3}; echo time (TE) of 3.8 ms; and repetition time (TR) of 2,300 ms. We collected dMRI using the following parameters of TE = 84 ms, TR = 9,000 ms, flip angle = 90°, reconstructed matrix size of 106 × 106, voxel size of 1.98 × 1.98 mm with slice thickness of 2 mm and a bandwidth of 1,745 Hz/Px. The data were recorded with 62 optimal nonlinear diffusion gradient directions at *b* = 1,500 s/mm^{2}. Approximately one nondiffusion weighted image (*b* = 0) per 10 diffusion-weighted images was acquired. We collected two dMRI datasets: one with anterior-to-posterior phase encoding direction and the other acquired in the opposite direction.

##### Tractography.

We used the structural connectivity between 90 AAL regions in the dMRI dataset described above. First, we used the linear registration tool from the FSL toolbox (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki; Functional Magnetic Resonance Imaging of the Brain [FMRIB], Oxford, UK) (48) to coregister the echo-planar imaging (EPI) image to the T1-weighted structural image. In turn, the T1-weighted image was coregistered to the T1 template of ICBM152 in Montreal Neurological Institute (MNI) space (49). We concatenated and inversed the resulting transformations, which were further applied to warp the AAL template (45) from MNI space to the EPI native space; ensured discrete labeling values were preserved by interpolation using nearest-neighbor method. The brain parcellations were thus conducted in each individual’s native space. The structural connectivity maps were generated for each participant using the dMRI data acquired. We combined the two datasets acquired with different phase encoding to optimize signal in difficult regions. Constructing these structural connectivity maps or structural brain networks involved a three-step process: 1) Regions of the whole-brain network were defined using the AAL template as used in the functional MRI data. 2) Probabilistic tractography was used to estimate the connections between nodes in the whole-brain network (i.e., edges). 3) Data were averaged across participants.

The FSL diffusion toolbox (Fdt) in FSL was used to carry out the various processing stages of the dMRI dataset. The default parameters of this imaging preprocessing pipeline were used for all participants. We then estimated the local probability distribution of fiber direction at each voxel (50). The probtrackx tool in Fdt was used to provide automatic estimation of crossing fibers within each voxel, which significantly improves the tracking sensitivity of nondominant fiber populations in the human brain (51).

We defined the connectivity probability from a given seed voxel *i* to another voxel *j* by the proportion of fibers passing through voxel *i* that reach voxel *j* when using a sampling of 5,000 streamlines per voxel (51). We then extended from the voxel level to the regional level, i.e., to a parcel in the AAL consisting of *n* voxels, i.e., 5,000 × *n* fibers were sampled. This allowed us to compute the connectivity probability *n* to region *p* as the number of sampled fibers in region *n* that connect the two regions divided by 5,000 × *v*, where *v* is the number of voxels in a given AAL region *i*.

We computed the connectivity probability from a given brain region to each of the other 89 regions within the AAL. It is important to note that the dependence of tractography on the seeding location, i.e., the probability from voxel *i* to voxel *j* is not necessarily equivalent to that from *j* to *i*. Still, these two probabilities are highly correlated across the brain for all participants (the least Pearson *r* = 0.70, *P* < 10^{−50}). We defined the unidirectional connectivity probability *n* and *p* by averaging these two connectivity probabilities since directionality of connections cannot be determined based on dMRI. This unidirectional connectivity can be thought of as a measure of the structural connectivity between the two areas, with

In order to consider the serotonin source region, we included in our model the connectivity between the raphe nucleus region and the rest of the brain. In this way, the neuronal activity across the whole brain excites through the existing fibers the raphe neurons, and this information is used in the neurotransmitter system (see below). In brief, we used the estimate of the raphe nucleus region from the *Harvard Ascending Arousal Network Atlas* (Athinoula A. Martinos Center for Biomedical Imaging; https://www.nmr.mgh.harvard.edu/resources/aan-atlas) (52). We then used *LeadDBS Mapper* (https://www.lead-dbs.org/) within Matlab to estimate the whole-brain tractography from this mask of the raphe nucleus to the rest of the brain (53, 54) using the Structural group connectome based on 32 Adult Diffusion HCP subjects GQI (55, 56). This allowed us to estimate the projections from the raphe nucleus to each AAL region, represented as a vector of 90 values.

It is important to note that the AAL parcellation has been shown to be less than optimal, in terms of being the least homogeneous parcellation scheme compared to a number of other parcellation schemes (57). Still, it is not clear whether the methodology used for this comparison is particularly meaningful. A recent paper by Eickhoff et al. (58) reviewed the literature on the topographic organization of the brain and concluded that there is no current consensus about what is the right spatial parcellation scheme. It is also important to note that, in contrast to these two papers, which are mostly concerned with the spatial organization of the brain, here we focus on the spatiotemporal global dynamics. As such, the AAL parcellation would seem a good choice given that AAL yields excellent significant results in the whole-brain literature in general (20, 35, 59), and, crucially, the relative low number of parcels in the AAL makes it highly suitable for our very extensive computational demands.

#### Functional dynamics.

The functional data are described in detail in a previously published study (13), but here we briefly summarize the participants, study setting, and acquisition protocol.

##### Ethics.

The study was approved by a National Health Service research ethics committee. Written informed consent was obtained from all participants prior to participation.

##### Participants in psilocybin study.

fMRI data from nine healthy subjects were included in this study (seven men; age, 32 ± 8.9 SD y of age). Study inclusion criteria were: at least 21 y of age, no personal or immediate family history of a major psychiatric disorder, substance dependence, cardiovascular disease, and no history of a significant adverse response to a hallucinogenic drug. All of the subjects had used psilocybin at least once before, but not within 6 wk of the study.

All subjects underwent two 12-min eyes-closed resting-state fMRI scans over separate sessions, at least 7 d apart. In each session, subjects were injected i.v. with either psilocybin (2 mg dissolved in 10 mL of saline, 60-s i.v. injection) or a placebo (10 mL of saline, 60-s i.v. injection) in a counterbalanced design. The injections were given manually by a medical doctor 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.

##### Anatomical scans.

Imaging was performed on a 3-T GE HDx system. These were 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. TR/TE, 7.9/3.0 ms; inversion time, 450 ms; flip angle, 20°.

##### BOLD fMRI data acquisition.

Two BOLD-weighted fMRI data were acquired using a gradient echo planer imaging sequence, TR/TE = 2,000/35 ms, field of view = 220 mm, 64 × 64 acquisition matrix, parallel acceleration factor = 2, 90° flip angle. Thirty-five oblique axial slices were acquired in an interleaved fashion, each 3.4 mm thick with zero slice gap (3.4-mm isotropic voxels). The precise length of each of the two BOLD scans was 7:20 min.

##### Preprocessing.

The fMRI data were preprocessed using MELODIC 3.14 (60), part of FSL (FMRIB’s Software Library; https://fsl.fmrib.ox.ac.uk/fsl/fslwiki) with standard parameters and without discarding any ICA components. For each participant and for each brain state (i.e., placebo and psilocybin), we used FSL tools to extract and average the BOLD time courses from all voxels within each region defined in the AAL atlas (considering only the 90 cortical and subcortical noncerebellar brain regions) (45).

##### Extracting PMSs.

Extracting the PMSs for the neuroimaging data relies on conducting the LEiDA analysis (31), summarized in Fig. 2. Similar to the procedure described in Deco et al. (29), we compute a phase coherence matrix to capture the amount of interregional BOLD signal synchrony at each time point, for all subject and conditions (placebo and psilocybin conditions). The phase coherence between each pair of brain regions is given by the following:

where the BOLD phases,

After computing the leading eigenvector of the phase coherence matrix *k*-means clustering algorithm with values of *k* from 2 to 20. Clustering solutions output a *k* number of cluster centroids in the shape of N×1 V_{C} vectors, which represent recurrent PMSs. The outer product of each cluster centroid vector *k* = 3 according many criteria including Silhouette analysis and the minimal *P* value for significant differences between probabilities between conditions.

Upon identifying PMSs, we computed the probability of occurrence and lifetime of each metastable substate in each condition. The probability of occurrence (or fractional occupancy) is given by the ratio of the number of epochs assigned to a given cluster centroid

Furthermore, to capture the trajectories of FC dynamics in a directional manner, we computed the switching matrix. This matrix contains the probabilities of transitioning from a given FC state (rows) to any of the other FC states (columns). Differences in probabilities of transition and probabilities of occurrence were statistically assessed between conditions using a permutation-based paired *t* test. This nonparametric test uses permutations of group labels to estimate the null distribution (computed independently for each experimental condition). To compare populations, we applied a *t* test to each of 1,000 permutations and used a significance threshold of α = 0.05.

#### Neurotransmitter receptor density.

A previous study has carefully described the methods used to obtain the 5-HT_{2A} receptor density distribution (11), but here we briefly summarize the main methods.

##### Ethics.

The study was approved by the Ethics Committee of Copenhagen and Frederiksberg, Denmark. Written informed consent was obtained from all participants prior to participation.

##### Participants.

The participants were healthy male and female controls from the freely available Cimbi database (63). The present data analysis was restricted to include individuals aged between 18 and 45 y. Acquisition of PET data and structural MRI scans were taken from 210 individual participants, yielding a total of 232 PET scans; of which 189 participants had only one scan, 20 participants had two scans, and a single had three scans.

##### PET and structural MRI.

Full details on PET and structural MRI acquisition parameters can be found in the original study (63) and in abbreviated form in ref. 9.

##### Extracting receptor density maps.

Similar to Deco et al. (9), we extracted the average receptor density for each individual AAL region *5-HT*_{1A}, *5-HT*_{1B}, *5-HT*_{2A}, and *5-HT*_{4} as well as *5-HTT* using standard FSL tools on the freely available receptor density maps in MNI space.

#### Whole-brain model neuronal system.

##### Whole-brain neuronal model (anatomy-plus-neuronal activity).

First, we modeled the system without any coupling between systems, i.e., a pure neurodynamical system (9). For this, we simulated the spontaneous brain activity at the level of the whole brain in a network model where each node represents a brain region and the links between nodes represent white matter connections. As proposed by Deco et al. (25), a dynamic mean field (DMF) model is used to simulate the activity in each brain region. In brief, based on the original reduction of Wong and Wang (64), this DMF model uses a reduced set of dynamical equations describing the activity of coupled excitatory (*E*) and inhibitory (*I*) pools of neurons to describe the activity of large ensembles of interconnected excitatory and inhibitory spiking neurons. The inhibitory currents, _{A} receptors, while excitatory synaptic currents, *n*, consists of reciprocally connected excitatory and inhibitory pools of neurons, whereas coupling between two areas *n* and *p* occurs only at the excitatory-to-excitatory level where it is scaled by the underlying anatomical connectivity *Materials and Methods*, *Structural Connectivity*).

The following system of coupled differential equations expresses the DMF model at the whole-brain level:

As can be seen, for each inhibitory (*I*) or excitatory (*E*) pool of neurons in each brain area *n*, the vector *E* and *I* pools is converted by the neuronal response functions, ^{−1} and ^{−1} determine the slope of H. Above the threshold currents of **6** and **7**, the uncorrelated standard Gaussian noise,

Emulating the resting-state condition, we used parameters in the DMF model based on Wong and Wang (64) such that each isolated node exhibited the typical noisy spontaneous activity with low firing rate *n,* we adjusted the inhibition weight,

The whole-brain network model used parcellated structural and functional MRI data from 90 cortical and subcortical brain regions (9), where each brain region *n* receives excitatory input from all structurally connected regions *p* into its excitatory pool, weighted by the underlying anatomical connectivity matrix, *Materials and Methods*, *Structural Connectivity*). Importantly, we used a global coupling factor *G* to equally scale all interarea E-to-E connections. Finding the optimal working point of the system simply requires finding the optimal global scaling factor, such that the simulated activity by the model is maximally fitting the empirical resting-state activity of participants in the placebo conditions.

The simulations were run with a time step of 1 ms for a range of *G* between 0 and 2.5 (with increments of 0.025). In order to emulate the empirical resting-state scans from nine participants, for each value of *G*, we ran 200 simulations of 435 s each.

Please note that to transform the simulated mean field activity from our DMF model into a BOLD signal, the generalized hemodynamic model of Stephan et al. (70) was used. This involves computing the BOLD signal in each brain area, *n*, arising from the firing rate of the excitatory pools

These biophysical variables are described by these equations:

where τ is a time constant, ρ is the resting oxygen extraction fraction, and the resistance of the veins is represented by *α*.

For each area, *n*, the BOLD signal,

We used all biophysical parameters stated by Stephan et al. (70) and concentrated on the most functionally relevant frequency range for resting-state activity, i.e., both simulated and empirical and BOLD signals were bandpass filtered between 0.1 and 0.01 Hz (71⇓⇓–74).

#### Mutually coupled whole-brain neuronal and neurotransmitter system.

Finally, in order to fully couple the whole-brain neuronal and neurotransmitter systems, we have to include an explicit description of the neurotransmitter dynamical system and the mutual coupling. We do this by modeling the dynamics of the neurotransmitter system through simulating the release-and-reuptake dynamics. These dynamics are then in turn coupled with the neuronal activity through the firing rate activity of the raphe brain region, source of the serotonin neurotransmitter.

Specifically, the concentration of serotonin is modeled through the Michaelis–Menten equations (26⇓–28):

where *n*, and where the corresponding density of a serotonin receptor **13** (i.e., brain raphe coupling). The **14**–**16**) is working in the nonlinear central sigmoidal part. For the second term on the right-hand side of Eq. **13**, ^{−1} is the maximum reuptake rate and

The reverse coupling, i.e., the effect of the neurotransmitter system on the neuronal system, is described in Eqs. **14**–**16**. Specifically, the neuronal activity is modeled as a dynamical system in Eqs. **14** and **15** (as a coupled variation of Eqs. **2** and **3**) by generating an extra current on the excitatory pyramidal and GABAergic inhibitory neurons in the following way:

where

where **14** and **15** describe the coupling of the neuronal to the neurotransmitter system. In this way, both neuronal and neurotransmitter dynamical system are explicitly expressed and mutually coupled.

#### Empirical fitting of mutually coupled whole-brain model to neuroimaging data.

We used the following measurements to measure the empirical fit of the mutually coupled whole-brain system.

##### Comparing empirical and simulated PMS space measurements.

We computed the symmetrized KLD between the simulated and empirical corresponding probabilities of the metastable substates, i.e., the probabilities of the extracted empirical centers after clusterization:

where *i*.

##### Comparing empirical and simulated transition probabilities between metastable substates.

We computed the entropy rate S of a Markov chain, with N states and transition matrix P. The rate entropy S is given by the following:

where

The stationary probability of state i is given by

The stationary distribution is the eigenvector of the transpose of the transition matrix with associated eigenvalue equal to 1. A Markov model with a lot of transitions will have a large rate entropy, while low entropy will be found in a Markov model with minimal transitions. For each transition matrix, we obtained the stationary distribution and, then, computed the entropy rate. Comparing the two transition probability matrices is defined by the absolute value of the difference between both respective Markov entropies.

### Data Availability.

The code to run the analysis and multimodal neuroimaging data from the experiment are available on GitHub (https://github.com/decolab/pnas-neuromod).

## Acknowledgments

G.D. is supported by the Spanish Research Project PSI2016-75688-P (Agencia Estatal de Investigación/Fondo Europeo de Desarrollo Regional, European Union); by the European Union’s Horizon 2020 Research and Innovation Programme under Grant Agreements 720270 (Human Brain Project [HBP] SGA1) and 785907 (HBP SGA2); and by the Catalan Agencia de Gestión de Ayudas Universitarias Programme 2017 SGR 1545. M.L.K. is supported by the European Research Council Consolidator Grant CAREGIVING (615539); Center for Music in the Brain, funded by the Danish National Research Foundation (DNRF117); and Centre for Eudaimonia and Human Flourishing, funded by the Pettit and Carlsberg Foundations. J. Cabral is supported by the Portuguese Foundation for Science and Technology (CEECIND/03325/2017), Portugal.

## Footnotes

- ↵
^{1}To whom correspondence may be addressed. Email: morten.kringelbach{at}psych.ox.ac.uk, nklmpg{at}stanford.edu, or gustavo.deco{at}upf.edu.

Author contributions: M.L.K. and G.D. designed research; M.L.K., J. Cruzat, J. Cabral, P.C.W., and G.D. performed research; M.L.K., G.M.K., R.C.-H., and G.D. contributed new reagents/analytic tools; M.L.K., J. Cruzat, J. Cabral, N.K.L., and G.D. analyzed data; and M.L.K., G.M.K., R.C.-H., P.C.W., N.K.L., and G.D. wrote the paper.

Reviewers: R.Q.Q., Leicester University; and O.S., Indiana University.

The authors declare no competing interest.

Database deposition: The code to run the analysis and the multimodal neuroimaging data (dMRI, fMRI, and PET) from the experiment are available on GitHub (https://github.com/decolab/pnas-neuromod).

See online for related content such as Commentaries.

- Copyright © 2020 the Author(s). Published by PNAS.

This open access article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

## References

- 1.↵
- 2.↵
- 3.↵
- G. Deco,
- J. Cruzat,
- M. L. Kringelbach

- 4.↵
- 5.↵
- 6.↵
- 7.↵
- 8.↵
- 9.↵
- 10.↵
- D. E. Nichols

- 11.↵
- V. Beliveau et al

- 12.↵
- L. D. Lord et al

- 13.↵
- R. L. Carhart-Harris et al

- 14.↵
- 15.↵
- 16.↵
- 17.↵
- 18.↵
- R. L. Carhart-Harris et al

- 19.↵
- R. L. Carhart-Harris et al

- 20.↵
- 21.↵
- 22.↵
- 23.↵
- C. J. Honey,
- R. Kötter,
- M. Breakspear,
- O. Sporns

- 24.↵
- G. Deco et al

- 25.↵
- G. Deco et al

- 26.↵
- 27.↵
- 28.↵
- L. Michaelis,
- M. L. Menten

- 29.↵
- G. Deco et al

- 30.↵
- 31.↵
- 32.↵
- 33.↵
- 34.↵
- 35.↵
- 36.↵
- 37.↵
- 38.↵
- A. Thiele,
- M. A. Bellgrove

- 39.↵
- A. Sajedin,
- M. B. Menhaj,
- A. H. Vahabie,
- S. Panzeri,
- H. Esteky

- 40.↵
- J. van Kempen,
- S. Panzeri,
- A. Thiele

- 41.↵
- R. L. Carhart-Harris

- 42.↵
- 43.↵
- 44.↵
- 45.↵
- 46.↵
- 47.↵
- 48.↵
- 49.↵
- 50.↵
- 51.↵
- 52.↵
- 53.↵
- 54.↵
- A. Horn et al

- 55.↵
- 56.↵
- 57.↵
- 58.↵
- S. B. Eickhoff,
- R. T. Constable,
- B. T. T. Yeo

- 59.↵
- 60.↵
- 61.↵
- 62.↵
- M. G. Preti,
- T. A. Bolton,
- D. V. Ville

- 63.↵
- G. M. Knudsen et al

- 64.↵
- K. F. Wong,
- X. J. Wang

- 65.↵
- 66.↵
- 67.
- 68.↵
- W. R. Softky,
- C. Koch

- 69.↵
- M. N. Shadlen,
- W. T. Newsome

- 70.↵
- 71.↵
- 72.↵
- 73.↵
- R. L. Buckner et al

- 74.↵
- S. Achard,
- R. Salvador,
- B. Whitcher,
- J. Suckling,
- E. Bullmore

## Citation Manager Formats

## Article Classifications

- Biological Sciences
- Neuroscience

## See related content: