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.

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. 4A). 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. 4B). 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. 4C) 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. 5B). 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. 5C) 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 mm3; 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/mm2. 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 P n p from region 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 P n p between regions 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 C n p = C p n . The computation of the regional connectivity probability was implemented using in-house Perl scripts. For both phase encoding directions, we constructed 90 × 90 symmetric weighted networks, based on the AAL90 parcellation, and normalized by the number of voxels in each AAL region. This represents the structural connectivity network organization of the brain. 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: dFC ( n , p , t ) = cos ( θ ( n , t ) − θ ( p , t ) ) , [1] where the BOLD phases, θ ( n , t ) , at a region n is estimated using the Hilbert transform for each BOLD regional time course. The Hilbert transform expresses a given signal x in polar coordinates as x ( t ) = A ( t ) ∗ cos ( θ ( t ) ) . Using the cosine function, two regions x and p with temporarily aligned BOLD signals (i.e., with similar angles) at a given TR will have a phase coherence value dFC ( n , p , t ) close to 1 [since cos ( 0 ° ) = 1 ]. On the other hand, time points where the BOLD signals are orthogonal (for instance, one increasing at 45° and the other decreasing at 45°) will have dFC ( n , p , t ) close to 0 [since cos ( 90 ° ) = 0 ]. The resulting dFC ( t ) for each participant in each condition is a 3D matrix with size N×N×T, where N = 90 is the number of brain regions and T is the total number of time points (here different for each subject and each condition). It is important to note that the phase coherence matrix is undirected and, as such, for each time t, the N×N dFC ( t ) matrix is symmetric across the diagonal. In order to characterize the evolution of the dFC matrix over time, we reduced the dimensionality of the problem by focusing on the evolution of the leading eigenvectors of the dFC matrices. The leading eigenvector, V 1 ( t ) , is a N×1 vector that captures the dominant connectivity pattern of the dFC ( t ) at each time t, given in matrix format by its outer product V 1. V 1 T . This approach substantially reduces the dimensionality of the data when compared to more traditional approaches considering all of the values in the connectivity matrix (43, 61, 62). Further details about LEiDA can be found in the work of Cabral et al. (31). After computing the leading eigenvector of the phase coherence matrix dFC ( t ) for each TR, we identified PMSs, i.e., recurrent FC patterns in the data, first described by Deco et al. (29). We detected a discrete number of FC patterns by clustering the leading eigenvectors V 1 ( t ) from the collapsed pre/post-psilocybin fMRI data (two conditions) including all participants. We ran the 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 V c . V c T is a N×N matrix representing the dominant connectivity pattern and the elements of V C weight the contribution of each brain area to the community structure established by the signs of the corresponding vector elements. We used HCP Workbench to render the cluster centroid vectors V C onto a cortical surface facilitate to visualization and interpretation of FC states. Here, we found that the optimal number of clusters was 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 V C divided by the total number of epochs (TRs) in each experimental condition. The lifetime is defined as the amount of time spent in each metastable substate (cluster centers). The probabilities and lifetimes were computed for each participant, in each experimental condition, and across the whole range of clustering solutions explored. 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, I ( I ) , are mediated by GABA A receptors, while excitatory synaptic currents, I ( E ) , are mediated by NMDA receptors. In terms of connectivity, each brain region, 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 C n p (Materials and Methods, Structural Connectivity). The following system of coupled differential equations expresses the DMF model at the whole-brain level: I n ( E ) = W E I 0 + w + J N M D A S n ( E ) + G J N M D A ∑ p C n p S p ( E ) − J n S n ( I ) , [2] I n ( I ) = W I I 0 + J N M D A S n ( E ) − S n ( I ) , [3] r n ( E ) = H ( E ) ( I n ( E ) ) = g E ( I n ( E ) − I t h r ( E ) ) 1 − exp ( − d E g E ( I n ( E ) − I t h r ( E ) ) ) , [4] r n ( I ) = H ( I ) ( I n ( I ) ) = g I ( I n ( I ) − I t h r ( I ) ) 1 − exp ( − d I g I ( I n ( I ) − I t h r ( I ) ) ) , [5] d S n ( E ) ( t ) d t = − S n ( E ) τ N M D A + ( 1 − S n ( E ) ) γ r n ( E ) + σ υ n ( t ) , [6] d S n ( I ) ( t ) d t = − S n ( I ) τ G A B A + r n ( I ) + σ υ n ( t ) . [7] As can be seen, for each inhibitory (I) or excitatory (E) pool of neurons in each brain area n, the vector I n ( E , I ) represents the total input current (in nanoamperes), the firing rate is denoted by the vector r n ( E , I ) (in hertz), while the synaptic gating is denoted by the vector S n ( E , I ) . The total input currents received by the E and I pools is converted by the neuronal response functions, H ( E , I ) , into firing rates, r n ( E , I ) , using the input–output function of Abbott and Chance (65), where the gain factors g E = 310 nC−1 and g I = 615 nC−1 determine the slope of H. Above the threshold currents of I t h r ( E ) = 0.403 nA, and I t h r ( I ) = 0.288 nA, the firing rates increase linearly with the input currents. The shape of the curvature of H around I t h r is given by the constants d E = 0.16 and d I = 0.087 . GABA receptors with τ G A B A = 0.01 s controls the average synaptic gating in inhibitory pools. NMDA receptors with a decay time constant τ N M D A = 0.1 s and γ = 0.641 / 1,000 (the factor 1,000 is for expressing everything in milliseconds) control the synaptic gating variable of excitatory pools, S n ( E ) . All excitatory synaptic couplings is weighted by J N M D A = 0.15 nA, while the weight of recurrent excitation is w + = 1.4 . The overall effective external input is I 0 = 0.382 nA with W E = 1 and W I = 0.7 . Note that in Eqs. 6 and 7, the uncorrelated standard Gaussian noise, υ n , has an amplitude of σ = 0.01 nA. 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 ( r ( E ) ∼ 3 Hz ) observed in electrophysiology experiments (66⇓⇓–69). Furthermore, similar to Deco et al. (25), for each node n, we adjusted the inhibition weight, J n , such that the firing rate of the excitatory pools r n ( E ) remained clamped at 3 Hz—even when receiving excitatory input from connected areas. The algorithm for achieving feedback inhibition control is described by Deco et al. (25), where it was shown that using this leads to a better prediction of the resting-state FC and to a more realistic evoked activity. 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, C n p , obtained from dMRI (9) (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 r n ( E ) , such that an increase in the firing rate causes an increase in a vasodilatory signal, s n , subject to autoregulatory feedback. In proportion to this vasodilatory signal, the blood inflow f n induces changes in deoxyhemoglobin content q n and blood volume v n . These biophysical variables are described by these equations: d s n / d t = 0.5 r n ( E ) + 3 − k s n − γ ( f n − 1 ) , [8] d f n / d t = s n , [9] τ d v n / d t = f n − v n α − 1 , [10] τ d q n / d t = f n ( 1 − ρ ) f n − 1 / ρ − q n v n α − 1 / v n , [11] 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, B n , is a static nonlinear function of deoxyhemoglobin content, q n , and volume, v n , that comprises a volume-weighted sum of extravascular and intravascular signals: B n = V 0 [ ( k 1 ( 1 − q n ) + k 2 ( 1 − q n v n ) + k 3 ( 1 − v n ) ] . [12] 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): d [ s n ] d t = α C B R r n ( E ) − V max [ s n ] ( K m + [ s n ] ) , [13] where [ s n ] is a vector, denoting the dynamics of the neurotransmitter concentration level in a brain region n, and where the corresponding density of a serotonin receptor R n is measured with PET. This is defining the interaction between the neuronal and neurotransmitter system through the first term on the right-hand side in Eq. 13 (i.e., brain raphe coupling). The C B R is outer product of the fiber density connectivity vector between the whole brain and the raphe nucleus, obtained through dMRI probabilistic tractography, and α = 5 is a factor normalizing the activity such that the current generated by the neurotransmitter (Eqs. 14–16) is working in the nonlinear central sigmoidal part. For the second term on the right-hand side of Eq. 13, V max = 1,300 nM⋅s−1 is the maximum reuptake rate and K m = 170 is the substrate concentration where the uptake proceeds at one-half of the maximum rate (26⇓–28). 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: I n ( E ) = W E I 0 + w + J N M D A S n ( E ) + G J N M D A ∑ p C n p S p ( E ) − J n S n ( I ) + W E S R n M n , [14] I n ( I ) = W I I 0 + J N M D A S n ( E ) − S n ( I ) + W I S R n M n , [15] where W E S and W I S are the excitatory and inhibitory feedback coupling parameters from the neurotransmitter system to the neuronal activity. The density of a serotonin receptor R n (as measured with PET) weighs the serotonin modulated current vector M n (current for each brain region), which is given by the sigmoid-like function commonly used in pharmacology and current equation (26): τ s d M n d t = − M n + J ( 1 + e − β ( l o g 10 [ s n ] + 1 ) ) , [16] where β = 10 , J = 0.1 , and τ s = 120 ms (26). The last terms on the right-hand side of Eqs. 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: K L ( P e m p , P s i m ) = 0.5 ( ∑ i P e m p ( i ) ln ( P e m p ( i ) P s i m ( i ) ) + ∑ i P s i m ( i ) ln ( P s i m ( i ) P e m p ( i ) ) ) , [17] where P e m p ( i ) and P s i m ( i ) are the empirical and simulated probabilities on the same empirical extracted metastable substates 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: S = S 1 + S 2 + … + S N , [18] where S i = − p ( i ) ∑ j = 1 N P ( i , j ) log P ( i , j ) . [19] The stationary probability of state i is given by p ( i ) . For long realizations of the Markov chain, the probabilities of each state converge to the stationary distribution p, which is the solution of the following equation: P T p = p . [20] 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 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.