A whole-brain model of the neural entropy increase elicited by psychedelic drugs

This whole-brain simulation reproduces earlier studies where increased entropy was found under the influence of psychedelics. The changes weren't uniform, with larger effects in the optical area. The changes were not well-explained by looking at where the most serotonin (2a) receptors were but mapped closely to the topological (spatial) properties of how the brain is connected.

Authors

Fernando Rosas
Enzo Tagliazucchi

Published

April 17, 2023

Scientific Reports

individual Study

Links

Read Paper DOI Google Scholar

Abstract

Psychedelic drugs, including lysergic acid diethylamide (LSD) and other agonists of the serotonin 2A receptor (5HT2A-R), induce drastic changes in subjective experience, and provide a unique opportunity to study the neurobiological basis of consciousness. One of the most notable neurophysiological signatures of psychedelics, increased entropy in spontaneous neural activity, is thought to be of relevance to the psychedelic experience, mediating both acute alterations in consciousness and long-term effects. However, no clear mechanistic explanation for this entropy increase has been put forward so far. We sought to do this here by building upon a recent whole-brain model of serotonergic neuromodulation, to study the entropic effects of 5HT2A-R activation. Our results reproduce the overall entropy increase observed in previous experiments in vivo, providing the first model-based explanation for this phenomenon. We also found that entropy changes were not uniform across the brain: entropy increased in all regions, but the larger effect were localised in visuo-occipital regions. Interestingly, at the whole-brain level, this reconfiguration was not well explained by 5HT2A-R density, but related closely to the topological properties of the brain’s anatomical connectivity. These results help us understand the mechanisms underlying the psychedelic state and, more generally, the pharmacological modulation of whole-brain activity.

Unlocked with Blossom Pro

Research Summary of 'A whole-brain model of the neural entropy increase elicited by psychedelic drugs'

Introduction

Psychedelic drugs acting at the serotonin 2A receptor (5HT2A-R), such as LSD, DMT and psilocybin, produce profound alterations of perception, cognition and selfhood. Earlier empirical work has identified two prominent neurophysiological signatures of the acute psychedelic state: suppression of alpha-band spectral power and an increase in the information-theoretic diversity of neural signals, commonly expressed as an increase in entropy. These entropy increases have been proposed to relate to both the acute alterations of consciousness and some longer-term psychological effects, and lie at the heart of the Entropic Brain Hypothesis (EBH), which links richness of subjective experience to the diversity of ongoing neural activity. Despite converging experimental evidence for elevated neural entropy under psychedelics, a clear mechanistic, whole-brain explanation for how 5HT2A-R activation produces these entropy changes has been lacking. Herzog and colleagues set out to provide a mechanistic account by extending a Dynamic Mean-Field (DMF) whole-brain model that incorporates neuromodulation by 5HT2A-R. The study aimed to test whether modulating neuronal response gain according to the empirical topography of 5HT2A-R could reproduce the experimentally observed increase in neural entropy, to map the regional patterning of entropy change, and to identify the roles of receptor density and structural connectivity (and their topological properties) in explaining those changes. By simulating resting-state population firing rates with and without 5HT2A-R activation and analysing regional differential entropy, the investigators sought a model-based explanation for the entropic effects attributed to serotonergic psychedelics.

Methods

The researchers used the Dynamic Mean-Field (DMF) model developed by Deco et al., which represents each brain region as coupled excitatory and inhibitory neural masses. Excitatory populations are coupled across regions via a structural connectivity matrix C, while local dynamics are governed by non-linear frequency–current (F–I) transfer functions. Neuromodulation was implemented as a multiplicative modulation of the F–I curve (a response-gain change) that scales with the regional density of the receptor of interest, d_rec_n. Key free parameters reported were the global coupling G and the neuromodulatory gain s_NM, set to 2.4 and 0.025 respectively following the authors’ fitting procedure. The model was integrated numerically using an Euler–Maruyama scheme and local feedback inhibition was tuned to maintain average firing rates near 3 Hz. Whole-brain simulations used the Automated Anatomical Labelling (AAL) parcellation (90 ROIs), a structural connectivity (SC) matrix averaged across 16 healthy subjects obtained from diffusion MRI tractography, and a 5HT2A-R density map derived from publicly available PET data. For each condition—placebo (PLA) and with 5HT2A-R activation—the DMF generated 90 excitatory firing-rate time series. The investigators treated 5HT2A-R activation as analogous to the pharmacological state elicited by serotonergic psychedelics. Differential entropy of each region's firing-rate time series was estimated by fitting a gamma distribution to the empirical firing-rate samples and computing the analytic closed-form expression for the gamma differential entropy from the fitted shape and scale parameters; goodness-of-fit by Kolmogorov–Smirnov statistics was reported as satisfactory. To explain spatial variation in entropy change, several linear mixed-effects models were fitted. The authors implemented an optimisation routine to partition regions into three groups (labelled S1, S2 and R1) according to whether entropy change was best predicted by connectivity strength (quadratic fit for strength-dependent regions) or by receptor density (linear fit for receptor-dependent regions). The final three-way model included fixed effects for connectivity strength and receptor density across the three groups and was evaluated by R2. To probe which topological features of the connectome drive the effects, the team generated three types of null-network surrogates that preserved increasingly strict attributes of the empirical connectome: (1) RAND — preserving overall density and strength distribution in a loose sense, (2) degree-preserving randomisation (DPR), and (3) strength-preserving randomisation (DSPR). For each surrogate network the DMF model was run with and without neuromodulation; each surrogate was simulated 120 times and results were averaged. Additional local graph measures (betweenness, eigenvector centrality, closeness, communicability, PageRank, subgraph centrality) were computed to test alternative predictors for regional entropy change.

Results

Simulating the DMF model with empirical 5HT2A-R topography produced a statistically significant increase in whole-brain differential entropy when comparing the 5HT2A-R condition to placebo. Mean global entropy rose from h_PLA = 2.15 nat to h_5HT2A = 2.25 nat (Wilcoxon signed-rank test p < 10^-6), with a large effect size (Cohen's d = 0.98 ± 0.05). Entropy increased across all 90 regions, but the magnitude of change was spatially heterogeneous: visuo‑occipital regions showed the largest increases, and regions within the primary visual and Default Mode (DMN) resting-state networks, as well as some frontoparietal executive network areas, exhibited pronounced entropy gains. At the whole-brain level, regional 5HT2A-R density was a poor predictor of the local entropy change (R2 = 0.0053 ± 0.0014). By contrast, connectivity strength (defined as the sum of weighted links for a region in the DTI-derived connectome) showed a moderate predictive relationship with local entropy changes (R2 = 0.45 ± 0.14). Visual inspection suggested two distinct tendencies in the strength–entropy relationship, motivating the optimisation-based partitioning of regions into three groups. This procedure yielded S1 (27 regions, mainly occipital, parietal and cingulate), S2 (42 regions), and R1 (11 regions, mainly temporal and frontal). For the S1 subset the strength–entropy association was very strong (S1: R2 = 0.93 ± 0.0032), whereas the R1 group’s entropy changes were strongly predicted by receptor density (R1: R2 = 0.93 ± 0.0085). Regions in S1 and S2 showed no significant relationship with receptor density. Combining connectivity strength, 5HT2A-R density and the three-way partition into a linear mixed-effects model accounted for 94.5 ± 0.002% of the variance in regional entropy change, indicating that these variables together provide a highly accurate statistical prediction of the modelled entropic effects. Analyses using other single-node connectivity measures (betweenness, eigenvector centrality, etc.) did not match the explanatory power of strength. Null-network experiments indicated that randomised and degree-preserving surrogate networks failed to reproduce the empirical pattern of entropy changes, while strength-preserving surrogates partially reproduced it. This pattern implies that the empirical distribution of node strengths is necessary but not sufficient to explain the full spatial reconfiguration; higher-order topological features beyond degree and strength distribution also contribute. Finally, the simulation findings aligned qualitatively with in vivo neuroimaging reports: occipital, parietal and cingulate entropy increases reported experimentally were mirrored in the model, and DMN regions showed increased entropy consistent with reduced DMN integrity observed under psychedelics.

Discussion

Herzog and colleagues interpret their findings as providing a mechanistic, whole-brain account that supports the Entropic Brain Hypothesis: modelling 5HT2A-R neuromodulation as a regionally varying gain modulation of neuronal response reproduced an average increase in neural entropy and a spatially heterogeneous reconfiguration across the cortex. Rather than a uniform amplification, 5HT2A-R activation produced larger entropy increases in highly connected sensory and associative regions—particularly occipital cortices and primary visual network nodes—suggesting a link between structural connectivity of specific regions and enhanced perceptual ‘bandwidth’ in the psychedelic state. The authors position their results relative to empirical studies that have reported occipital, parietal and cingulate entropy increases under LSD and psilocybin, noting that many such regions fell into the strength-dependent S1 group in their model. They further highlight that increased entropy in DMN regions in the simulations is consistent with experimental reports of diminished DMN integrity and with phenomenology such as ego‑dissolution. From a mechanistic viewpoint, the study suggests that psychedelics act not merely locally by receptor binding but by amplifying collective, emergent dynamics shaped by the topology of anatomical connectivity; connector hubs and highly anatomically connected regions may therefore have outsized influence on global dynamical regularity. Several limitations are acknowledged. The structural connectome derived from diffusion MRI is incomplete and averaged across subjects, and the DMF model treats each region as internally homogeneous, omitting fine-grained local circuitry that may underlie phenomena like geometric visual hallucinations. The simulations considered only 5HT2A-R neuromodulation, whereas classic psychedelics interact with additional receptors (for example dopaminergic receptors in the case of LSD). These simplifications may preclude reproducing other empirical signatures of psychedelics such as alpha suppression or directed functional-connectivity changes. The authors also note that their analyses focused on univariate regional statistics and did not characterise multivariate interactions or information flow across regions. For future work the paper suggests several directions advocated by the investigators: applying forward models to map simulated firing rates to fMRI or M/EEG observables for more direct comparisons with empirical datasets; exploring non-linearities and dose–response relationships in the model; incorporating individual-level connectome and receptor maps to build personalised predictive models; and extending the analysis to multivariate information-theoretic measures (including approaches related to Integrated Information Theory) to better capture network-level phenomena thought to underlie complex subjective effects. The authors conclude that their results constitute a first mechanistic explanation for the neural entropy increase induced by 5HT2A-R activation and emphasise the need to move from a unidimensional to a multidimensional characterisation of consciousness in psychedelic research.

View full paper sections

SECTION

short-and long-term effects of the psychedelic experience, including certain aspects of the reported subjective experienceand subsequent personality changes. Interestingly, the opposite effects have been reported for states of loss of consciousness, where a strong decrease in brain entropy has been repeatedly observed. This seems to be a core feature of loss of consciousness, generalising across states such as deep sleep, general anaesthesia, and disorders of consciousness. Together, these facts are yielding converging evidence that entropy and related measures offer simple and powerful indices of consciousness. Based on these findings, Carhart-Harris and colleagues have put forward the Entropic Brain Hypothesis (EBH) an entropy-based theory of conscious states. The EBH proposes the simple, yet powerful idea that the variety of states of consciousness can be indexed by the entropy of key parameters of brain activity; or, in other words, that the richness of subjective experience is directly related to the richness of on-going neural activity, where richness can be defined most simply as diversity. Note that even when different states can be indexed by their corresponding diversity, we are not claiming that one state is more conscious than other. Therefore, investigating the neurobiological origins of such changes in brain activity is a key step in the study of altered states of consciousness. Multiple experiments in humans and animal models have established that the mind-altering effects of psychedelics depend on agonism specifically at the serotonin 2A receptor (5HT2A-R). A recent simulation study involving whole-brain computational modelling confirmed that the topographic distribution of 5HT2A-R in the human brain is critical to reproduce the functional connectivity dynamics (FCD) of human fMRI data recorded under the acute effects of LSD. Here we build upon this model to characterise the interplay between the entropy of brain signals, the distribution of 5HT2A receptors, and structural connectivity properties of the brain, with the overarching goal of explaining the sources of entropic effects, and thus altered consciousness, elicited by psychedelic drugs.

RESULTS

We simulated whole-brain activity using the Dynamic Mean-Field (DMF) model by Deco et al., using parameter values fit to reproduce the dynamics of fMRI recordings in humans during wakeful rest, as well as under the acute effects of LSD (see Supplementary Figure). The model consists of interacting pools of excitatory and inhibitory neural populations, coupled via long-range excitatory connections informed by the anatomical connectivity of the brain. The DMF model combines a theoretical model of neural and synaptic dynamics with two empirical sources of information about the human brain: the human connectome, i.e. DTI-estimated connectivity between the 90 regions of the AALatlas; and average 5HT2A-R expression across the human brain obtained with Positron Emission Tomography (PET) scans. Each simulation of the DMF model generates 90 time series of excitatory and inhibitory firing rates, one of them for each region of the AAL atlas (Fig.). These firing rates have a non-linear dependency on the local inputs, determined by a frequency-current (F-I) curve. Given that 5HT2A-R seems to modulate the sensitivity of neurons rather than driving them towards excitation or inhibition, we followed Deco et al.and modelled 5HT2A-R activation as a response-gain modulationof this F-I curve dependent on the receptor density at a given region (Fig.). We simulated the model in two conditions, with and without 5HT2A-R activation, which, in an analogy with neuroscientific terminology, we refer to as placebo (PLA) and 5HT2A-R conditions, respectively. In this way, we obtain 90 time series of excitatory firing rates for each condition, which we keep for further analysis. Finally, using these time series, we estimate Shannon's differential entropy for each region under both conditions, yielding a topographical distribution of entropy values (Fig.) that constitutes the main subject of analysis in this study. Further details about model specification and entropy estimation, as well as other methodological caveats, are presented in the Methods section at the end of this article.

HT2A-R ACTIVATION CAUSES A HETEROGENEOUS, LINEAR INCREASE IN THE ENTROPY OF SIMULATED BRAIN ACTIVITY.

In our first analysis, we used the DMF model to test the main prediction of the EBH: that 5HT2A-R activation causes an increase in the overall entropy of neural signals (Fig.). In line with previous experimental findings with psychedelic drugs, the model shows a significant increase in the brain's entropy h due to 5HT2A-R activation, with an average entropy of h PLA = 2.15 nat in the placebo condition, and h 5HT2A = 2.25 nat in the 5HT2A-R condition (dashed vertical lines in Fig., Wilcoxon signed-rank test p < 10 -6 , Cohen's d = 0.98 ± 0.05). A closer look at the distribution of entropy changes, however, reveals a more nuanced picture, entropy increases for all regions, but the effect is larger for visuo-occipital regions, compared to the rest of the regions (Fig., Supp. Fig.). This suggests that, according to the model, 5HT2A-R agonism might trigger a complex reconfiguration of the topographically specific distribution of entropy, and not a mere homogeneous overall increase. To investigate this reconfiguration, we analysed the local changes in entropy by plotting the entropy of the nth region in both PLA and 5HT2A conditions, denoted by h 5HT2A n and h PLA n respectively (Fig.). Our results show that 5HT2A-R activation affects local entropy in a non-uniform and linear manner, especially in regions with high baseline resting state entropy. In particular, in such regions the effect of 5HT2A-R activation show a different tendency than the rest of the regions, such that the local entropy in the 5HT2A-R condition could not always be determined by the region's baseline entropy. This finding hints towards a more general theme: that local dynamical properties alone are not able to explain the local changes in activity induced by 5HT2A-R agonists like psychedelic drugs. We will explore this phenomenon in depth in the following sections. Next, we studied the effect of 5HT2A-R agonism on local entropy by considering the relative change scores, Based on recent in vivo experiments with serotonergic psychedelics we expect to find localized entropy increases on occipital, cingulate and parietal regionsas well as several changes on regions belonging to Resting State Networks (RSNs). Accordingly, we split brain regions by standard anatomical and functional groupings, and found that the effect were localised specially in occipital areas, together with some cingulate and parietal regions. (Fig., Supp. Fig.). Additionally, regions participating in the primary visual and Default Mode (DMN) RSNs-including occipital and cingulate areas, respectively-showed a marked tendency to increase their entropy, as well as regions participating in the Frontoparietal Executive Network (FPN). Further discussion about the relation between these results and recent in vivo experiments with psychedelic drugs is presented in the Discussion section. To explain the effects of 5HT2A-R activation, the first natural step is to factor in the density of 5HT2A-R at each region. Somewhat surprisingly, at a whole-brain level, receptor density was a very poor predictor of the entropy change due to 5HT2A-R activation (Fig., R 2 = 0.0053 ± 0.0014 ). In contrast, we found that the connectivity strength (based on the DTI human connectome), defined as the sum of all the weighted links connecting a given region, exhibited a moderate correlation with local entropy changes (Fig., R 2 = 0.45 ± 0.14). By visual inspection, however, it is clear that the relationship between DTI-informed connectivity strength and h n is complex and can be split in two separate tendencies. To investigate this phenomenon, we implemented a simple optimisation algorithm to find the set of regions with the strongest relationship between strength and h n . This yielded a set of regions with a higher strength-h n correlation (S 1 , R 2 = 0.93 ± 0.0032 , as opposed to R 2 = 0.45 ± 0.14 for the whole-brain), and induced a partition of all brain regions into three groups: those Changes in entropy were overall independent from receptor density, although (B) they were well predicted by the connectivity strength of each region. We split into strength (blue and gray), and receptor dependent groups (red). The S 1 and S 2 groups showed no significant relationship with receptor density, while the R 1 group were highly correlated with it. (C) Topographical localisation of the three groups, following the same colour code. S 1 were mainly located in occipital, parietal and cingulate regions, while the R 1 ones were in temporal and frontal ones. highly dependent on strength (S 1 , 27 regions, blue in Fig.), those less dependent on strength (S 2 , 42 regions , gray), and those highly dependent on receptors (R 1 , 11 regions, red). Building up on this partition, we studied the effect of receptor density on those regions where entropy change could not be predicted from connectivity strength-the R 1 group. We found a strong relationship between receptor density and entropy changes for the R 1 group, which resulted in positive correlation ( R 2 = 0.93 ± 0.0085 ). On the contrary, receptor density does not predict entropy changes for areas in the S 1 and S 2 . This shows a complementary role of density and strength in S 1 , S 2 , and R 1 regions: entropy changes in R 1 regions strongly depend on the receptor density, but not connectivity strength; while changes in S 1 , S 2 regions depend on connectivity strength, but not receptor density. Overall, to assess the predictive power of connectivity strength and the receptor density on local entropy changes, we built a linear mixed model for h n using connectivity strength, 5HT2A-R density, and the afore- mentioned three-way separation of brain regions as the predictor variables. Together these variables explain 94.5 ± 0.002% of the variance of h n , confirming that they provide an accurate model for predicting the entropic effects of 5HT2A-R activation. This suggests that psychedelic drugs and other 5HT2A-R agonists do not have a simple, localised effect on brain activity, but instead amplify the fundamentally collective, emergent properties of the brain as a complex system of interacting elements. The specific connectivity strength distribution explains relative changes in entropy. As a final analysis, we set out to investigate exactly which topological properties of the brain's structural connectivity explain the observed changes in entropy. With this exercise, we aim to answer two questions: whether any property simpler than connectivity strength can explain the results; and whether any property more complicated than connectivity strength is needed to do so. To this end, we ran further simulations of the DMF model using suitable null network models of the human connectome (Fig.), with the 5HT2A-R density map held fixed. In particular, we used three null models designed to test increasingly strict null hypotheses that preserved different network attributes of the original connectome: (1) the overall density and strength (RAND, Fig.), (2) the degree distribution (degree-preserving randomisation [DPR], Fig.); and (3) the strength distribution (strength-preserving randomisation, [DSPR], Fig.). We simulated the DMF model in these surrogate networks, with and without neuromodulation, and computed the resulting entropy changes h n . Our first result is that random and degree-preserving surrogate networks are unable to reproduce the entropy changes observed in the human connectome: when compared against the h n values obtained from the unper- turbed network, neither of them produce values close to the original (Fig.). This result asserts the findings in the previous section, showing that indeed node strength plays an important role in shaping the global pattern of entropy change associated with the action of psychedelics and other 5HT2A-R agonists. Simpler topological features, like the degree distribution, are not enough to explain such changes. Perhaps more interestingly, the connectivity strength-preserving surrogate networks partially reproduced the entropy changes, suggesting that higher-order properties of the connectome are necessary (Fig.). Furthermore, we repeated the analyses in the previous section trying to predict h n using other local connectivity measures (e.g. betweenness, eigenvector centrality), and none of them produced a model as accurate as strength (Supp. Fig.). Together, these results show that, once the receptor distribution is fixed, the network strength distribution is both necessary but not sufficient to explain the entropic effects of 5HT2A-R activation. As previously suggested, other topological network features are known to mediate transitions of consciousness in other contexts, and investigating which network properties explain high-order dynamical signaturesof psychedelics remains an exciting avenue of future work.

DISCUSSION

In this study, we investigated the brain entropy changes induced by serotonergic psychedelics by simulating whole-brain resting state activity with and without 5HT2A-R activation. In contrast to empirical studies, which usually only have access to coarse-grained fMRI or M/EEG data, our approach allows us to study firing rates of brain regions and hence to directly assess the effect of 5HT2A-R agonism at the level of neural population activity. In agreement with the Entropic Brain Hypothesis, the model shows a significant increase of global brain entropy, with a inhomogeneous and linear reconfiguration. The diversity of effects stresses the importance of extending the scope of the EBH from a simple level-based approach towards a multi-dimensional perspective, which might better characterise the richness of 5HT2A-R activation effect both on the brain and consciousness. We propose a possible mechanism in the sense that the average increase in neural entropy induced by serotonergic psychedelic drugs can be explained (in statistical sense) in terms of quantities that provide a neurobiological interpretation about how this phenomena could occur in the brain. These results have important consequences for our understanding of consciousness, neurodynamics, and the psychedelic state. In the following, we highlight some of those implications, and propose possibilities for future work. 5HT2A-R-induced entropy changes are regionally heterogeneous. Our main result in this study is a mechanistic explanation of the main prediction of the EBH that serotonergic psychedelics increase the entropy of brain activity. However, one of the main takeaways is that this overall increase is tied to a spatially heterogeneous reconfiguration, rather than a globally consistent increase, in entropy. Our simulation results show that 5HT2A-R activation triggers an entropy increase in sensory areas, as observed in occipital cortices, as well as the primary visual RSN (Supp. Fig.). This agrees with the increased perceptual 'bandwidth' that characterises the psychedelic state, as higher entropy might be related to richer perceptual experience in a given moment of consciousness-potentially related to reduced gating. As a matter of fact, all the primary visual RSN regions are part of the S 1 group, suggesting that the effect of psychedelics on perceptual experience might be directly related to the connectivity of those regions. The localisation of the entropy increases may relate to domain-specific changes in consciousness, which could be interpreted as consistent with a recent dimensions of consciousness characterisation of the psychedelic experience. Comparison to in vivo experiments with psychedelic drugs. Throughout the paper, we have focused on one particular signature of 5HT2A-R agonists on the brain-a global increase in average entropy. But much more is known about the effect of psychedelics on the brain, and studying these more nuanced effects is key to understanding the rich phenomenology of the psychedelic state. In this section we deepen this connection by providing a more detailed comparison between the behaviour of the model and experimental results with psychedelic drugs. Our findings are in agreement with earlier studies where the effect of psychedelic drugs on the topographical distribution of entropy-related measures was correlated with subjective effects. For example, Schartner et al.studied the effect of LSD, psilocybin, and ketamine on the entropy rate of binarized MEG signals, and found localised increases in entropy rate on occipital and parietal regions. Similarly, Lebedev et al.analysed the sample entropy of fMRI recordings of subjects under the effect of LSD, finding localized increases on frontoparietal, medial occipital, posterior and dorsal cingulate regions. Many of those regions showed a consistent increase of their entropy with 5HT2A-R agonism in our study (Fig.). Moreover, many of those regions actually belong to the S 1 group (c.f. Fig.), which suggests that their entropy increase in experimental data might be directly related to the connectivity in those regions. Together, these findings support the conclusion that the DMF model, once optimised, can reproduce not only functional connectivity, but also some of the most salient localised entropy increases found on in vivo human studies up to date. Also, on a more fundamental level, our findings suggest that psychedelics disrupt the functional organisation of the brain with an especially focal and pronounced action on highly anatomically connected brain regions. In fact, experimental and modelling evidence points towards a key role of connector hubs on psychedelic and psychopatological states. Perturbing the activity in these specific regions (via 5HT2A-R agonism) could have particularly profound implications for the regularity of brain activity and the quality of conscious experience. However, sub cortical regions may also play an important role on the modulation of brain activity, conscious experience, and the psychedelic state. On a separate line of inquiry, there is strong evidence associating the DMN to high-level cognitive functions, including the sense of space, time and recognition of (self) boundaries. Disruptions to the DMN have been linked to fundamental changes in experience, such as ego dissolution. Our simulation results show that 5HT2A-R activation increases the entropy of all DMN regions, which is consistent with the reported decrease in the DMN network integrityinduced by psychedelic drugs. In contrast, low-level motor functions such as motor regulation remain largely unaffected during the psychedelic state, which is consistent with the modest entropy changes observed in the lingual and superior motor areas.

CURRENT LIMITATIONS AND FUTURE RESEARCH.

The approach employed in this work presents certain limitations related to several aspects of the simulation and analysis. Acknowledging and understanding these limitations can help us extend and improve our approach, while introducing new questions in the field of psychedelic computational neuroscience. To our knowledge, the DMF model is the only whole-brain model that implements neuromodulation and is capable of reproducing neuroimaging data in the placebo and psychedelic states. Nonetheless, it makes some important simplifications that are worth discussing. At the network level, the DTI-based connectome used here is known to be incomplete, thus improvements could be made to the model parameters of brain connectivity. At the dynamical level, the DMF model models neuronal populations as perfectly homogeneous within a given region, and it is known that finer-grained local structure of certain brain regions is likely to be key to explaining certain subjective effects of the psychedelic state (e.g. lattice structure in the visual cortex and geometric visual hallucinations). Additionally, the version of the DMF model used here only considers 5HT2A-R agonism, while classic serotonergic psychedelic drugs also have high binding affinities for other receptors (e.g. in the case of LSD, the D1 and D2 dopamine receptors). These simplifications do not prevent the model from reproducing statistical features of brain signals under the placebo and 5HT2A-R conditions, but could result in an inability to reproduce finer aspects of the dynamics of the whole-brain activity in these conditions. Extending the model to reproduce other dynamical signatures of psychedelics (like alpha suppressionor reduced directed functional connectivity) constitutes a natural extension of this work, with the recent example of Kringelbach et al.where an extension of the DMF can reproduce dynamical features of the brain in placebo and in the psychedelic state Another potentially fruitful line of future work involves making more detailed comparisons with in vivo psychedelic neuroimaging data, and, potentially, subjective experience reports. For example, one natural option would be to use forward models of fMRIor M/EEGto bridge between the firing rates produced by the DMF model and other data modalities, to produce simulated data that is more directly comparable with available empirical data. Another exciting possibility is to explore model parameters to examine potential non-linearities and their implications for different relevant aspects of brain function. For example, there are reasons to believe that the dose-response relationship is non-linear for psychedelics and that over a certain threshold dosage (level of 5HT2A-R stimulation) new subjective and global brain function properties can appear. Most interestingly, a potentially very useful extension of this work would be to include individual subjectlevel connectome and receptor data to build personalised models of response to psychedelic drugs. This would enable a much more comprehensive modelling framework, capable of correlating structural brain features with subjective experience reports. Such a framework could potentially make individualised predictions of the action of serotonergic psychedelics on specific individuals, aiding patient stratification and treatment customisation. Finally, it is worth noting that all our analyses here are based on the univariate statistics of individual brain regions, not including any correlation or information flow between them. However, it is known that some highlevel subjective effects of psychedelics (such as complex imageryand ego dissolution) are related to network, as opposed to single region, dynamics. Therefore, building a richer statistical description of the brain's dynamics using recent information-theoretic tools (such as multivariate extensions of mutual informationor Integrated Information Theory, IIT) remains an exciting open problem. In fact, recent attempts of unifying IIT and the EBH in a single framework to understand the effects of psychedelic drugsare helping to bridge the gap between the univariate EBH and the multivariate IIT. Final remarks. In this paper, we have provided the first mechanistic explanation of the neural entropy increase elicited by psychedelic drugs, using a whole-brain dynamical model with 5HT2A-R neuromodulation. Furthermore, we built a simple model able to predict a region's relative change in entropy from its local 5HT2A-R density and topological properties, showing that, somewhat paradoxically, at a whole-brain level receptor density is a poor predictor of 5HT2A-R activation effect. Key to developing this predictive model was a three-way partition of brain regions according to the dependency they exhibited with connectivity strength, suggesting a differentiated action mechanism of 5HT2A-R agonists that depends on the local topology of brain regions. In summary, our results suggest that the local changes in entropy, as well as the global entropy increase, induced by 5HT2A-R activation can be explained from a region-specific interplay between structural connectivity and receptor density. Finally, controlled experiments with null network models confirm that receptor density and connectivity strength are necessary, but not sufficient, to explain the entropic effects of 5HT2A-R activation. The spatially heterogeneous, complex nature of the observed effects of 5HT2A-R activation opens a challenging problem for understanding the clinical and scientific relevance of psychedelic drugs and their entropic effect. Furthermore, it stresses the necessity of moving beyond the current unidimensional approach to consciousness to a multi-dimensional one, that better captures the phenomenological and neurodynamical richness of psychedelic state.

METHODS

Dynamic mean-field model with 5HT2A-R neuromodulation. The main computational tool used in this study is the Dynamic Mean-Field (DMF) model by Deco et al.which consists of a set of coupled differential equations modelling the average activity of one or more interacting brain regions. In this model, each brain region n is modelled as two reciprocally coupled neuronal masses, one excitatory and one inhibitory, and the excitatory and inhibitory synaptic currents I (E) and I (I) are mediated by NMDA and GABA A receptors respectively. Different brain regions are coupled via their excitatory populations only, and the structural connectivity is given by the matrix C. The full model, including the neuromodulatory effect described below, is given by Above, for each excitatory (E) and inhibitory (I) neural mass, the quantities n , and S (E,I) n represent its total input current (nA), firing rate (Hz) and synaptic gating variable, respectively. The function F(•) is the transfer function (or F-I curve), representing the non-linear relationship between the input current and the output firing rate of a neural population. Finally, J FIC n is the local feedback inhibitory control of region n, which is optimized to keep its average firing rate at approximately 3 Hz (Supp. Fig 1), and v n is uncorrelated Gaussian noise injected to region n. The interested reader is referred to the original publication for further details. Key to this study is the recent extension of this model including neuromodulatory effects. In the equations above, g NM n is a neuromodulatory scaling factor modulating the F-I curve of all brain regions in the model, which affects region n depending on its density of the receptor of interest, d rec n . Including neuromodulation, the free parameters of the model are G, the global coupling parameter, and s NM , the neuromodulatory gain, that are set to 2.4 and 0.025, following our fitting procedure (Supp. Fig.). To simulate the expression of 5-HT2A receptor on the excitatory and inhibitory populations, neuromodulation was applied to both populations. The model was simulated using a standard Euler-Maruyama integration method, using the parameter values shown in Table. Table. Dynamic mean field (DMF) model parameters. Parcellation, connectome, and 5HT2A receptor maps. In addition to the DMF equations above, to simulate whole-brain activity we need three more ingredients: a suitable parcellation of the cortex into welldefined regions, a structural connectivity matrix between those regions, and the density map of a receptor of interest-in our case, the serotonin 2A receptor. We clarify that all the empirical inputs of the model were produced elsewhere. For a detailed description of the procedure, please follow the corresponding references cited below. As a basis for our simulation, we used the Automated Anatomical Labelling (AAL) atlas, a sulci-based parcellation of brain volume registered in MNI space. The AAL parcellation specifies a partition of the brain into 90 Regions Of Interest (ROIs), that provides sufficient level of detail to obtain a picture of the topographical distribution of entropy, while being suitable for computationally intensive simulations. Structural connectivity data between the 90 AAL ROIs was obtained from 16 healthy subjects (5 female) using diffusion Magnetic Resonance Imaging (dMRI), registered to the MNI space and parcellated according to the AAL atlas on the subject's native space. Then, for each subject the histogram of fiber directions at each voxel was obtained, yielding an estimate of the number of fibers passing trough each voxel. The weight of the connection between regions i and j was defined as the proportion of fibers passing through region i that reach region j, which (since the directionality of connections cannot be determined using dMRI) yields a symmetric 90 × 90 Structural Connectivity (SC) matrix. Finally, the SC matrix of each subject was thresholded, pruning any connection lower than 0.1%, and SC matrices of all subjects were averaged to obtain a single SC matrix used in all simulations. This average matrix was kindly provided by et al.For a detailed description of these methods see or refer to the Supplementary Subsection 7. For the density map of 5HT2A receptors, we used the Positron Emission Tomography (PET) data made public by Beliveau et al., which can be combined with the AAL parcellation to obtain an estimate of receptor density in every AAL region. Further details can be found in the original publications. Together with the DMF equations, these three elements fully specify our whole-brain model and allow us to run simulations and obtain time series of excitatory firing rates, r (E) n , which we save for further analysis. Estimating the differential entropy of firing rates. The result of each simulation is a set of 90 time series representing the firing rate of each excitatory population ( r (E) n ), which we analysed to measure the entropy of the activity of each ROI. For a continuous random variable X with associated probability distribution p(x) and support X , its differential entropy is defined asEstimating differential entropy from data is in general a hard problem, stemming (for the most part) from the difficulties in estimating a probability density from samples. We found that the probability distributions of firing rates r (E) n were well approximated by a gamma distribution: goodness-of-fit, measured with a standard Kolmogorov-Smirnov statistic, was satisfactory and comparable across brain regions and conditions (see Supp. Fig.). This observation greatly simplifies the estimation of the differential entropy of each region, which is now reduced to estimating the shape and scale parameters of the gamma distribution, k and θ . Once these parameters are estimated (in our case by standard maximum likelihood estimation), the differential entropy of a gamma distribution can be computed analytically in closed form as where Ŵ(•) and ψ(•) are the standard gamma and digamma functions, respectively. Linear models for h n and three-way separation of AAL regions. To explain the changes in dif- ferential entropy across regions and conditions, we trained several linear mixed-effects modelsusing different sets of covariates designed to test specific hypotheses about the relation between entropy, 5HT2A-R density, and network topology. First, we implemented a simple algorithm to find the group of regions where the change in entropy depends on region strength. The model was fit by maximum likelihood using off-the-shelf software, and explanatory power measured using R 2 . To do so, we designed an optimisation procedure to find the optimal fit between h n and connectivity strength, using all occipital regions of the AAL90 atlas as a prior for a quadratic fit (i.e. the S 1 group). The rest of the regions were used for a linear fit between h n and receptor density. Regions were assigned to the S 1 group if the residuals of the quadratic fit to strength were lower than the residuals for the linear fit to receptor density. An optimization algorithm randomly partitioned the remaining regions in two groups until a maximum in the R 2 of a linear mixed-effects model with 3 grouping variables (i.e. S 1 , S 2 , and R 2 ) is obtained. The model used a quadratic fit for strength-dependent regions, and a linear one for receptor-dependent regions. This three-way separation of brain regions induces a model with 7 terms: a constant term and 2 fixed-effect terms for each group of regions: one for the connectivity strength and one for the receptor density. Null network models of the human connectome. Null network models of the human structural connectivity were used to evaluate the role of the local connectivity propertieson the local changes in entropy induced by 5HT2A-R activation. To this end, we applied three different randomisation schemes to the structural connectivity in order to produce suitable surrogate networks. These surrogates were designed to preserve different network attributes of the original connectome: (1) the overall density and strength (RAND), (2) the degree h(X) = - After randomisation, the DMF model was run with and without 5HT2A-R activation, and entropies estimated and analysed following the same procedure as in the rest of the article. Every surrogate model was run 120 times, and the results averaged across runs. In addition to the network surrogates, we computed several topological measures to include in the linear model reported in Fig.. These measures included betweenness, eigenvector centrality, closeness, communicability, PageRank index, and sub-graph centrality. All of these computations, including the surrogate networks and the topology measures, were performed using the Brain Connectivity Toolbox.