RETRACTED ARTICLE: A mechanistic model of the neural entropy increase elicited by psychedelic drugs

Psychedelic drugs, including lysergic acid diethylamide 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 some regions and decreased in others, suggesting a topographical reconfiguration mediated by 5HT2A-R activation. 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.


Supplementary Subsection 1. Heterogeneous entropy changes induced by 5HT2A-R activation on anatomical and functional groupings of brain regions.
We asked whether the heterogeneous changes of entropy induced by 5HT2A-R activation could be explained by grouping the AAL brain regions according to anatomical and functional criteria. Regarding anatomical criteria brain regions can be spatially split into 8 major non-overlapping anatomical groups: Frontal, Temporal, Parietal, Occipital, Limbic, Sensorimotor, Cingulate and Subcortical. 1 Regarding the functional grouping, we used the Resting State Networks, a functional grouping of brain regions based on the observed spatio-temporal patterns of BOLD signals during resting state activity. 2 These groups are Salience (Sal), Fronto-Parietal (FPN), Default Mode (DMN), Primary Visual (Vis), Extrastiate Cortex (EC), Auditory (Aud), Sensorimotor (SM), and Executive Control (EC). The first one, Sal, was obtained from Leeet al. 3 , the second and third, FPN and DMN, were obtained from Oliver et al. 4 , and the rest from Beckmann et al. 2 Brain regions can potentially belong to different functional groups.
The effect of 5HT2A-R activation on entropy is heterogeneous also at the level of anatomical groups (Supplementary Figure 1A) -i.e. within all the groups there were both regions with increased and decreased entropy after the 5HT2A-R activation. However, Occipital and Cingulate regions show a strong tendency to increase their entropy, in agreement with entropy increases in these regions observed in human experiments with serotonergic psychedelics. 5,6 Regarding the functional grouping (Supplementary Figure 1B), we also found that the effect of 5HT2A-R activation on regional entropy (Fig. 3B) is heterogeneous within groups, with the exception of Vis and FPN, where almost all the regions increased and decreased their entropy, respectively. Note also that with the exception of the angular gyri, all the DMN regions increase their entropy, which resonates with the observed reduction of the DMN integrity on humans during psychedelic experiences. 7 Figure 1. The effect of 5HT2A-R activation is not explained by anatomical or functional groupings. (A) 5HT2A-R activation heterogeneously change entropy in regions belonging to the same anatomical group. (B) Entropy also changed in a heterogeneous way when functional grouping is considered. Note the special case of FPN and Vis networks, for which almost all regions decreased and increased their entropy, respectively (note that functional groups can be overlapping). Circles are averages and error bars 1 s.d. computed from 1000 simulations.

Supplementary Subsection 2. Connectivity strength is the best predictor for entropy changes among local connectivity measures.
We control the role of local connectivity on the observed entropy changes (∆h n ) induced by 5HT2A-R activation using as predictors for ∆h n other local connectivity measures as: degree, 8 eigenvector centrality, 8 communicability, 9 page-rank, 8 sub-graph centrality, 8 and closeness centrality. 10 We confirmed that local connectivity strength is the best linear predictor of entropy changes among other centrality measures. Figure 2. Centrality measures as possible explanatory variables for ∆h n . Local connectivity strength is the best linear predictor of ∆h n among other centrality measures. Communicability and sub-graph centrality also show a good linear relationship with ∆h n , though with higher residual variance. Circles are averages and error bars 1 s.d. out of 1000 independent simulations.

Supplementary Subsection 3. Receptors controls.
We control the role of receptor density distribution on the observed entropy changes (∆h n ) induced by 5HT2A-R activation by using alternative receptor distributions. To test the effect of the heterogeneous 5HT2A-R distribution on ∆h n , we used an uniform receptor distribution, where the receptor density for each region corresponds to the average 5HT2A-R density among all regions. To evaluate the role of the specific 5HT2A-R distribution, a randomised version of the 5HT2A-R was also used, where for each simulation the receptor expression is randomised. Finally, to evaluate the effect of activating other serotonin receptors, we used the 5HT4-R distribution, a receptor thought to play no role on the psychedelic state.
On the one hand, the Uniform and Random receptor distributions show almost the same behaviour on average, with ∆h n showing a linear relationship with the node strength that changed slope at node strength value ∼ 0.35. The Random distribution shows more variability because the receptor distribution is randomised on each simulation. On the other hand, the 5HT4-R distribution changes entropy but in a subtler manner, yielding an average ∆h n ∼ 0.
On summary, the heterogeneous 5HT2A-R distribution critically impacts the heterogeneity of ∆h n values and its relationship with the local connectivity strength. In addition, using other empirically-based distribution yields no significant change in average brain entropy, highlighting the key relationship between 5HT2A-R density and connectivity on the entropic effects induced by serotonergic psychedelic drugs. Figure 3. Controls for receptor distribution vs ∆h n . ∆h n vs connectivity strength is plotted for 5HT2A-R, Uniform, 5HT4-R, and Random distribution of receptors (see text for details). The specific distribution of 5HT2A-R is critical both for the heterogeneity of ∆h n and for the average increase in local entropy. Circles are averages and error bars 1 s.d. out of 1000 independent simulations.

Supplementary Subsection 4. The probability distribution of excitatory firing rates can be well fitted by a Gamma distribution.
To check the goodness of fit (GOF) of Gamma distribution to the simulated excitatory firing rates of each region, we generated 10 6 simulation points for each region under both PLA and 5HT2A condition, then, a Gamma distribution was fitted and the Kolmogorov-Smirnov (K-S) distance between the firing rate distribution and 1000 random samples (same size) of the respective fitted Gamma distribution was computed. We summarise the Gamma GOF as the average K-S for each region under both conditions. To assess the significance of average GOF values, we generated an acceptance interval computing the K-S distance between 100 independent samples (same size) sampled from exactly the same Gamma distribution (re-sampled K-S). This procedure was repeated for different set of Gamma distribution parameters. If the average K-S distance of a given region falls below the maximum re-sampled K-S distance, we consider this region to be well fitted by the Gamma distribution.
Despite the decreased GOF under 5HT2A condition (compared to PLA) all regions fall within the acceptance interval under both conditions, enabling us to use the Gamma distribution parameters to estimate each region's Shannon's differential entropy.   For the sake of completeness, in the following we describe in detail the methods used by Ref. 11 to obtain the human connectome. We clarify that we did not produce these connectome data sets, instead, we used data kindly provided by Gustavo Deco.

Supplementary
The structural connectivity between the 90 AAL regions was obtained from averaging across 16 healthy young adults (5 females, mean SD age: 24.75 2.54). The linear registration tool from the FSL toolbox (www. fmrib.ox.ac.uk/fsl, FMRIB, Oxford) was used to coregister the EPI image to the T1-weighted structural image. The T1-weighted image was co-registered to the T1 template of ICBM152 in MNI space. The resulting transformations were concatenated and inversed and further applied to warp the AAL template from MNI space to the EPI native space, where interpolation using nearest-neighbor method ensured that the discrete labeling values were preserved. Thus the brain parcellations were conducted in each individual's native space. The structural connectivity (SC) maps were generated for each participant using the dMRI data acquired. The two data sets acquired have different phase encoding to optimize signal in difficult regions. The construction of these structural connectivity maps consisted of a three-step process. First, the regions of the whole-brain network were defined using the AAL template as used in the functional MRI data. Second, the connections between nodes in the whole-brain network (i.e., edges) were estimated using probabilistic tractography. Third, data was averaged across participants.
The FSL diffusion toolbox (Fdt) was used to carry out the various processing stages of the diffusion MRI data using the default parameters of this imaging pre-processing pipeline on all participants. Following this preprocessing, we estimated the local probability distribution of fiber direction at each voxel. The probtrackx tool in Fdt was used to provide automatic estimation of crossing fibers within each voxel.
The connectivity probability from a seed voxel i to another voxel j was defined by the proportion of fibers passing through voxel i that reach voxel j using a sampling of 5000 streamlines per voxel. This was extended from the voxel level to the region level, i.e., in an AAL parcel consisting of n voxels, 5000x n fibers were sampled. The connectivity probability P i j from region i to region j is calculated as the number of sampled fibers in region i that connect the two regions divided by 5000 x n, where n is the number of voxels in region i. The SC matrix was thresholded at 0.1%, i.e., five streamlines.
For each brain region, the connectivity probability to each of the other 89 regions within the AAL was calculated. Due to the dependence of tractography on the seeding location, the probability from i to j is not necessarily equivalent to that from j to i. However, these two probabilities are highly correlated across the brain for all participants (the least Pearson r = 0.70, p < 10 − 50). As directionality of connections cannot be determined based on diffusion MRI, the unidirectional connectivity probability P i j between regions i and j was defined by averaging these two connectivity probabilities. This unidirectional connectivity was considered as a measure of the structural connectivity between the two areas, with C i j = C ji . The regional connectivity probability was calculated using in-house Perl scripts. For both phase encoding directions, 90x90 symmetric weighted networks were constructed based on the AAL90 parcellation, and normalized by the number of voxels in each AAL region; thus representing the structural connectivity network organization of the brain.