# Simulated 5-HT2A receptor activation accounts for the high complexity of brain activity during psychedelic states

Hugo M Martin, Rodrigo Cofre, Alain Destexhe

► To cite this version:

Hugo M Martin, Rodrigo Cofre, Alain Destexhe. Simulated 5-HT2A receptor activation accounts for the high complexity of brain activity during psychedelic states. 2026. hal-05473680

HAL Id: hal-05473680

https://inria.hal.science/hal-05473680v1

Preprint submitted on 23 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Simulated 5-HT2A receptor activation accounts for the high complexity of brain activity during psychedelic states

Hugo M. Martin<fnref n="1" />, Rodrigo Cofre<fnref n="1" />,<fnref n="2" />, Alain Destexhe<fnref n="1" />*

## Abstract

Serotonergic psychedelics, such as LSD, psilocybin, and DMT, have strong effects on human brain activity, yet their mechanisms of action at the whole-brain level are only partially understood. Here, we present a biophysically-based mean-field model that integrates cellular and network-level details to simulate the effects of these compounds at different spatial scales. By incorporating the brain-wide distribution of 5-HT<sub>2A</sub> receptors, our model mechanistically links receptor activation to a reduction in leak membrane potassium conductance, consistent with electrophysiological data. Our simulations reveal that this microscopic perturbation leads to the emergence of a brain state characterized by asynchronous and irregular dynamics with increased firing rates, as well as significant alterations in spectral power. Specifically, we find a robust decrease in power within the delta, theta, and alpha frequency bands, a result consistent with empirical findings. This change in dynamics is accompanied by an increase in spontaneous complexity, as quantified by the Lempel-Ziv complexity index, as observed experimentally. Furthermore, our model accurately replicates experimental findings regarding the Perturbational Complexity Index (PCI), demonstrating that PCI does not increase significantly by psychedelic drug administration. This crucial dissociation, where spontaneous complexity and spectral power are increased while perturbational complexity is preserved, highlights the distinct neurophysiological substrates underlying different metrics in psychedelic states. Our multiscale model provides a robust, mechanistic framework for understanding how serotonergic psychedelics modulate global brain activity, offering new insights consistent with empirical neuroimaging and electrophysiological data.

## 1 Introduction

Over the past few years, research on psychedelics has gained significant momentum, driven predominantly by their therapeutic potential in treating mood disorders and addictions [1, 2, 3]. This growing interest emphasize the need to mechanistically understand the effects of psychedelics compounds on the brain and the need too integrate in biophysical models consistent findings from the experimental literature.

Recent advances in neuroimaging have provided novel insights into brain organization across multiple scales, enabling the development of computational models that integrate biophysical mechanisms with large-scale neural dynamics [4, 5, 6, 7, 8]. These whole-brain models provide a rich tool for identifying robust biomarkers of different states of consciousness, offering a principled framework to investigate how changes at the molecular and cellular scales translate into macroscopic patterns of brain activity [9]. In this context, modeling the action of serotonergic psychedelics at the whole-brain level provides a unique opportunity to explore the dynamical principles underlying consciousness alterations, allowing the integration of different neuroimaging

data such as functional and diffusion MRI and PET scans [10, 11, 12, 13, 14, 15].

The action of serotonergic psychedelics in the brain originates at the molecular and cellular levels, where receptor binding and intracellular signaling cascades drive alterations in neuronal excitability and network activity [16, 17, 18, 19]. The integration of molecular and cellular details into whole-brain models remains computationally prohibitive, emphasizing the need for theoretically grounded frameworks that link molecular scale interactions to whole brain neural dynamics.

In this work, we use a computational framework designed to bridge the gap between molecular and cellular processes and whole-brain dynamics, addressing the challenge of biologically plausible large-scale brain simulations. At the heart of our approach is a biophysically principled mean-field model that incorporates cellular and molecular mechanisms while reducing the system's dimensionality to enable computationally efficient large-scale simulations. By capturing emergent network properties across multiple spatial scales, our model enables realistic whole-brain simulations of serotonergic psychedelic action. Our hierarchical modeling approach is based in a previously published framework [9], that consists of four interdependent levels:

<footnote n="1">¹ Paris-Saclay University, CNRS, Paris-Saclay Institute of Neuroscience (NeuroPSI), 91400 Saclay, France</footnote>

<footnote n="2">² Université Côte d'Azur, INRIA CRONOS Team, Sophia Antipolis, France</footnote>

<footnote>*Corresponding author: alain destexhe@cnrs.fr</footnote>
Simulated 5-HT2A receptor activation accounts for the high complexity of brain activity during psychedelic states

(i) single-cell models that describe molecular interactions and intrinsic neuronal dynamics (microscale);
(ii) networks of spiking neurons that characterize local circuit activity (mesoscale); (iii) a mean-field model that approximates collective network behavior while maintaining biophysical accuracy (mesoscale, low-dimensional); and (iv) a whole-brain model leveraging mean-field dynamics to simulate large-scale neural activity constrained by empirical anatomical connectomes ( macroscale). To implement and validate our simulations, we employ The Virtual Brain (TVB) platform [20].

Specifically, here we present a bottom-up model of serotonergic psychedelic action that explicitly incorporates the effects of reduced leak potassium (K+) conductance. The proposed mean-field model integrates membrane conductances and synaptic receptor dynamics to simulate serotonergic modulation mediated by 5-HT₂A receptors, which are considered the principal targets through which serotonergic psychedelics elicit hallucinogenic effects in humans [21, 22]. Activation of these receptors has been shown to modulate leak K+ conductances, consistent with electrophysiological evidence [23, 24, 25, 26, 27, 28].

At the whole brain scale, psychedelic compounds consistently increase spontaneous signal diversity in the human brain, a neural marker associated with altered conscious content. Studies using electroencephalography (EEG) and magnetoencephalography (MEG) have shown increased neuronal complexity in a variety of psychedelic states. These effects are robustly captured by metrics such as Lempel–Ziv complexity (LZc) and are frequently accompanied by reductions in delta-, theta- and alpha-band oscillatory power [29, 30, 31, 32, 33, 34, 35].

Importantly, such increases in spontaneous complexity do not necessarily translate to changes in evoked cortical dynamics. Studies employing transcranial magnetic stimulation coupled with EEG (TMS–EEG) report that while psychedelics such as psilocybin and ketamine enhance spontaneous signal diversity, the Perturbational Complexity Index (PCI)—an index of the brain’s causal capacity to integrate and differentiate information—does not increase significantly [32, 31]. This dissociation suggests that spontaneous complexity may reflect the richness of conscious content, whereas evoked complexity may be more tightly linked to the underlying capacity for consciousness itself.

Incorporating receptor density distributions across the brain, our model predicts that serotoninergic psychedelics induce an activated brain state characterized by asynchronous-irregular dynamics and increased complexity relative to the waking resting

state, linking synaptic receptor activity to whole-brain network dynamics.

Consistent with empirical findings, we observe out model exhibit heightened brain complexity, as measured by LZc. Furthermore, our simulations of evoked activity show that, in agreement with experimental data, PCI increases, but not significantly under psychedelics. We also observe a robust decrease in spectral power within the $\delta$, $\theta$, and $\alpha$ bands under reduced K+ conductance, consistent with the literature on serotonergic psychedelics.

## 2 Results

This section introduce the framework that connects single-neuron excitability to large-scale activity. We then show how the decreased potassium leak conductance increases neural excitability and, in turn, drives measurable alterations in brain state, reflected in spontaneous complexity (LZc), evoked complexity (PCI), and spectral power.

### 2.1 Multiscale Modeling Framework from Single Neurons to Whole-Brain Dynamics

To investigate how cellular-level mechanisms shape whole-brain dynamics under serotonergic modulation, we implemented a multiscale modeling framework spanning single neurons, local circuits, and the macroscale connectome (Figure 2).

At the microscale, we employed the Adaptive Exponential integrate-and-fire (AdEx) model, balancing biological plausibility with computational tractability in large-scale simulations [36, 37]. The AdEx model captures a wide repertoire of cortical firing behaviors, including regular spiking and fast spiking. The model is governed by two coupled differential equations: one describing membrane potential dynamics (Eq.7) and another describing the adaptation current (Eq.8), which regulates neuronal excitability as a function of spiking history. Different parameterizations allow the representation of excitatory pyramidal neurons and inhibitory interneurons, fitted from in vitro electrophysiological recordings [38]. This versatility permits adaptation to multiple brain structures, including the thalamus, basal ganglia, and cerebellum [39, 40, 41].

At the mesoscale, we constructed networks of spiking neurons to represent local microcircuits or anatomically defined regions. In a canonical cortical column model, 10,000 neurons with an 80:20 ratio of pyramidal cells to interneurons which are randomly connected with probability $p = 0.05$, consistent with empirical estimates [42]. Synaptic transmission was
Simulated 5-HT2A receptor activation accounts for the high complexity of brain activity during psychedelic states

<figure><img src="image_2.png" /></figure>

**Figure 1:** 5-HT<sub>2A</sub> receptor activation depolarizes cortical neurons by reducing K<sup>+</sup> leak conductance. (A) Schematic illustrating the biophysical mechanism where 5-HT<sub>2A</sub> receptor activation closes potassium (K<sup>+</sup>) leak channels, reducing outward current and leading to depolarization. (B) Electrophysiological traces from regular-spiking (RS, top) and fast-spiking (FS, bottom) cortical neurons in control conditions (left) and following simulated 5-HT<sub>2A</sub> receptor activation via reduced K<sup>+</sup> leak conductance (right). For both cell types, reducing K<sup>+</sup> leak conductance from 8.21 to 5.37 nS depolarizes the resting membrane potential from -70 to -55 mV and increases firing rate in response to the same current injection. This models the excitatory effect of serotonergic modulation mediated by 5-HT<sub>2A</sub> receptors.

modeled as conductance-based via AMPA/NMDA (excitatory) and GABA<sub>A</sub> (inhibitory) receptors.

To efficiently bridge from high-dimensional spiking networks to population dynamics, we employed a mean-field approach derived from a second-order Master Equation [43]. This approach preserves key biophysical features—such as mean activity, conductance fluctuations, and adaptation—while enabling rapid simulation of large-scale systems. The core of this method is the neuronal transfer function, which relates firing rate output to presynaptic input. For conductance-based models without closed-form solutions, we use a semi-analytical procedure: numerically simulating neuron responses across input conditions and fitting them with analytical templates, validated for AdEx neurons [44, 38] and extended to more complex neuron types such as the quadratic integrate-and-fire or Hodgkin-Huxley formulations [45, 46, 41].

At the macroscale, the brain was modeled as a

network of mean-field units coupled according to empirically derived structural connectivity. The chosen parcellation defines the model's spatial resolution. We used the 68-region Desikan-Killiany atlas [47] for the results in the main article, but we replicate all our results using the cortically-reduced Lausanne 463 parcellation. Inter-regional coupling strengths and conduction delays were derived from diffusion MRI tractography, with delays computed from anatomical distances. Each regional unit receives both stochastic background input, modeled as Ornstein-Uhlenbeck process, and weighted, time-delayed inputs from other regions.

Figure 2 illustrates this hierarchical architecture. The top row shows the three scales: single-neuron AdEx models reproducing diverse firing patterns and modulations, including those induced by 5-HT<sub>2A</sub> receptor activation and altered K<sup>+</sup> leak conductance; local populations represented either as explicit excitatory-inhibitory spiking networks or as mean-field
Simulated 5-HT2A receptor activation accounts for the high complexity of brain activity during psychedelic states

models capturing their mean dynamics; and a whole-brain network linking mean-field units according to the anatomical connectivity. The bottom row shows the application of this framework to simulate both spontaneous and evoked conditions (additional external stimulus is applied to a single region) under different physiological regimes. From spontaneous activity, we extracted measures such as Lempel-Ziv complexity (LZc) and spectral features; from evoked activity, we computed the Perturbational Complexity Index (PCI).

## 2.2 5-HT<sub>2A</sub> receptors and K+ conductance

Activation of 5-HT<sub>2A</sub> receptors leads to the closure of K+ channels (see Figure 1), resulting in neuronal depolarization, increased excitability, and a higher firing rate (see Supplementary Fig. S1). Below we show how we implement this in the model.

**Leak Term Decomposition:** The dynamics of an isopotential neuron's membrane potential are governed by:

$$C \frac{dV}{dt} = -g_L(V - E_L) + F(V, t) \quad (1)$$

where:

* C: Membrane specific capacitance (F/cm²)

* V: Membrane potential (mV)

* g<sub>L</sub>: Leak conductance (nS)

* E<sub>L</sub>: Leak reversal potential (mV)

* F(V, t): Function encompassing synaptic inputs and active currents

To model the action of 5-HT<sub>2A</sub> receptors, we take into account experiments showing that the stimulation of these receptors blocks leak K+ conductances in all cortical neuron types [23, 24, 25, 26, 27, 28]. To do this, we express the leak conductance g<sub>L</sub> as a sum of Na+ and K+ leak conductances, g<sub>Na</sub> and g<sub>K</sub>, such as:

$$g_L(V - E_L) = g_{Na}(V - E_{Na}) + g_K(V - E_K) \quad (2)$$

where E<sub>Na</sub> and E<sub>K</sub> are their respective reversal potentials.

Expanding and matching terms in (2):

Conductance balance: $g_L = g_{Na} + g_K$ \quad (3)

Reversal potential: $E_L = \frac{g_{Na}E_{Na} + g_K E_K}{g_{Na} + g_K}$ \quad (4)

From (3) and (4), we derive:

$$g_{Na} = g_L - g_K \quad (5)$$

Solving (4) for g<sub>K</sub>:

$$g_L E_L = g_{Na} E_{Na} + g_K E_K \\ g_L E_L = (g_L - g_K) E_{Na} + g_K E_K \quad (\text{from (5)}) \\ g_K = g_L \frac{E_L - E_{Na}}{E_K - E_{Na}} \quad (6)$$

**Baseline State:** Given the following baseline values:

* g<sub>L</sub> = 10 nS, E<sub>L</sub> = -65 mV

* E<sub>Na</sub> = 50 mV, E<sub>K</sub> = -90 mV

Using (6), we compute the following:

$$g_K = 10 \frac{-65 - 50}{-90 - 50} \approx 8.21 \text{ nS}$$

Then from (5):

$$g_{Na} = 10 - 8.21 \approx 1.79 \text{ nS}$$

For g<sub>L</sub> = 10 nS and E<sub>L</sub> = -65 mV, we have g<sub>Na</sub> ≈ 1.79 nS and g<sub>K</sub> ≈ 8.21 nS.

**Depolarized State:** To model the effect of 5-HT<sub>2A</sub> receptors, we can reduce g<sub>K</sub> until the resting membrane potential reaches a desired value, here -55 mV (see Fig. 1). To compute the new g<sub>K</sub> associated with this target value, we fix g<sub>Na</sub> to its previously computed value (here 1.79 nS) and solve the previous system for g<sub>L</sub> and g<sub>K</sub>:

$$g_L = g_{Na} \frac{E_{Na} - E_K}{E_L - E_K} \quad \text{(rearranged from (4))} \\ = 1.79 \frac{50 - (-90)}{-55 - (-90)} = 7.16 \text{ nS}$$

$$g_K = g_L - g_{Na} \quad \text{(from (3))} \\ \approx 7.16 - 1.79 = 5.37 \text{ nS}$$

When 5-HT<sub>2A</sub> receptors are activated, g<sub>K</sub> decreases, leading to an increase in E<sub>L</sub> and V and, consequently, an increase in average firing rate (see Supplementary Fig. S1).

## 2.3 Lempel-Ziv Complexity as a Marker of Consciousness in Psychoactive States

A central hypothesis in consciousness research proposes that conscious experience arises from the brain's capacity to generate and integrate complex patterns of neural activity [48, 49]. This complexity can be quantified using measures such as Lempel-Ziv complexity (LZc), a non-parametric metric of algorithmic complexity that estimates the number of distinct subsequences in a time series, thereby indexing its information content and unpredictability [50].

In EEG and MEG studies, LZc has proven sensitive to state-dependent changes in neural dynamics,
<figure><img src="image_3.png" /><figcaption>**Figure 2: Multiscale Modeling Framework for Brain Dynamics.** Our modeling framework integrates neural activity across three scales: single neurons, local populations, and the whole brain. **Top row:** At the single neuron scale, we used the Adaptive Exponential (AdEx) integrate-and-fire model, because its capacity to reproduce diverse firing patterns such as regular spiking and fast spiking, influenced by factors like 5-HT<sub>2A</sub> receptor activation and K<sup>+</sup> leak currents. The population scale illustrates two complementary approaches: a detailed spiking neural network (left), depicting excitatory (blue) and inhibitory (red) populations with external inputs (V<sub>ext</sub>), and a mean-field model (right), showing the mean dynamics of excitatory (Exc) and inhibitory (Inh) populations and their interactions. At the whole brain scale, the model is built by connecting multiple mean-field units based on anatomical connectivity. **Bottom row:** The framework's utility in simulating different conditions (e.g., control and 5-HT<sub>2A</sub> activation) under both spontaneous and evoked simulation paradigms. Spontaneous activity can be analyzed for metrics such as Lempel-Ziv complexity (LZc) and spectral properties. Evoked simulations are used to analyze responses to external stimuli, leading to the quantification of the Perturbational Complexity Index (PCI). This hierarchical approach enables the investigation of how cellular-level mechanisms give rise to emergent macroscopic brain states.</figcaption></figure>

with particular relevance to serotonergic psychedelic states [51, 30, 32, 33]. Accumulating evidence indicates that these agents reliably increase LZc, reflecting heightened signal diversity and information richness during altered states of consciousness. Increases in LZc are consistently associated with subjective reports of intensified perceptual and cognitive experiences, supporting the “entropic brain” hypothesis [52, 53], which posits that the richness of conscious experience emerges from a dynamic balance between order and disorder in brain activity. Conversely, reduced-consciousness states such as deep non-REM sleep, general anesthesia, or

coma are characterized by decreased LZc, consistent with synchronized, low-complexity neural dynamics [54, 55, 56, 57, 58, 59, 60, 61, 62].

While LZc effectively captures global differences in neural dynamics, interpreting its changes under psychedelics requires a mechanistic bridge from molecular-level 5-HT<sub>2A</sub> receptor agonism to whole-brain effects. To address this, we developed a biophysically based whole-brain model that simulates cortical network dynamics under serotonergic modulation. We modeled the effects of reduced potassium (K<sup>+</sup>) conductance—consistent with 5-HT<sub>2A</sub>-mediated depolarization—on network complexity. Specifically,
Simulated 5-HT2A receptor activation accounts for the high complexity of brain activity during psychedelic states

<figure><img src="image_4.png" /></figure>

**Figure 3: Lempel-Ziv complexity (LZc) and Robustness Analysis.** Left) This figure shows the difference in Lempel-Ziv complexity (LZc) between the Baseline and reduced K+ leak conductance conditions. The box plot displays the median (red line), interquartile range (box), and whiskers extending to 1.5 times the interquartile range. The violin plot shows the data distribution, and individual data points correspond to whole brain simulations are overlaid. A significant difference, as indicated by the three asterisks (t-test, *p* < 0.001), is observed between the two conditions, with the reduced K+ leak conductance condition showing a higher LZc value. This is consistent with an increase in complexity reported in the literature under serotonergic psychedelics. Right) The heatmap demonstrates the robustness of the LZc results to changes in the target excitatory and inhibitory reversal potentials (baseline EL<sub>e</sub> = -63 and EL<sub>i</sub> = -65). The highlighted cell in green corresponds to the parameters used in the figure on the left. The color bar on the right shows the normalized LZc difference, where red indicates a positive difference (increase in LZc) and blue a negative difference. The asterisks indicate the level of statistical significance (t-tests, ***: *p* < 0.001, **: *p* < 0.01, *: *p* < 0.05, ns: *p* > 0.05). The prevalence of significant differences across a wide range of parameter values confirms that the observed increase in LZc under the reduced K+ leak conductance condition is a robust finding.

we set our model in two conditions, a baseline condition and a reduced K+ leak conductance described in Section 2.2 and run 15 seeds of spontaneous activity of each condition. Figure 3 (left) shows that K+ conductance reduction significantly increased simulated LZc relative to baseline, mirroring the elevated complexity reported in empirical EEG/MEG data under psychedelics. This effect was robust, as confirmed by a systematic sensitivity analysis (Figure 3, right) in which LZc differences persisted across a broad range of increased excitatory and inhibitory reversal potentials. The prevalence of significant increases across the parameter space demonstrates that the modeled LZc enhancement is not contingent on fine-tuned conditions but emerges as a general property of the network under 5-HT<sub>2A</sub>-like modulation. These results provide mechanistic support for the hypothesis that serotonergic psychedelics enhance cortical complexity, aligning both with empirical findings and theoretical accounts of an entropic reorganization of brain dynamics. To further probe the specificity of this effect, we examined the distribution of LZc across regions under various receptor maps (see Supplementary Figs. S2 and S3), as depicted in boxplots for baseline conditions alongside 5-HT<sub>2A</sub> receptor-mediated implementations, which applied region-

specific reductions in leak potassium conductance proportional to empirical 5-HT<sub>2A</sub> receptor density maps.

## 2.4 Perturbational Complexity Index in Simulated Baseline and Psychedelic States

The Perturbational Complexity Index (PCI) quantifies the brain's capacity to integrate and differentiate information in response to a controlled external perturbation [63]. It is computed from the algorithmic complexity of transcranial magnetic stimulation (TMS)-evoked electroencephalographic (EEG) responses and has been validated as a robust, theory-informed biomarker of the conscious level, independent of task demands or sensory processing. High PCI values are a hallmark of wakefulness, reflecting rich and distributed cortical interactions, whereas markedly reduced PCI values consistently characterize states of diminished consciousness, including non-rapid eye movement (NREM) sleep, general anaesthesia, and disorders of consciousness such as the vegetative state/unresponsive wakefulness syndrome (VS/UWS) [63, 64].

Recent empirical findings have revealed a notable dissociation between PCI and the subjective
Simulated 5-HT2A receptor activation accounts for the high complexity of brain activity during psychedelic states

<figure><img src="image_5.png" /></figure>

**Figure 4: Perturbational Complexity Index (PCI) Analysis and Robustness.** Left) The box plot compares the Perturbational Complexity Index (PCI) values between the Baseline and reduced K+ leak conductance conditions. The central red lines represent the median PCI values, while the boxes indicate the interquartile range (IQR). The whiskers extend to 1.5 times the IQR. The violin plots show the distribution of the simulations and individual points are overlaid. The results show no statistically significant difference ("ns", t-test, *p* > 0.05) in PCI between the two conditions, which is consistent with findings in the literature regarding the effects of psilocybin on PCI. Right) The heatmap demonstrates the robustness of the PCI results under changes in key model parameters. The x- and y-axes represent variations in the excitatory and inhibitory target reversal potentials (baseline *EL*<sub>e</sub> = -63 and *EL*<sub>i</sub> = -65). The highlighted cell in green corresponds to the parameters used in the figure on the left. The color scale represents the normalized PCI difference, with red indicating a positive difference and blue a negative difference. The abundance of "ns" (not significant) labels across the heatmap confirms that the lack of a significant change in PCI, showing the robustness of our findings across a wide range of these parameter values.

alterations of consciousness induced by serotonergic
psychedelics: despite marked changes in perception
and cognition, PCI values remain in the range typical
for normal wakefulness [32]. A similar effect has
been reported under low-dose ketamine, which pro-
duces a non-serotonergic psychoactive state[31]. This
suggests that while psychedelics profoundly modify
the content of consciousness, they do not diminish
its capacity, in contrast to unconscious states. Such
results highlight PCI's specificity to global integra-
tive capacity and motivate its combined use with
spontaneous signal complexity measures—such as
Lempel–Ziv complexity—to capture the full neuro-
physiological impact of psychedelics.

To mechanistically investigate this dissociation, we
used our model in two conditions: (i) baseline wake-
fulness, and (ii) a psychedelic state implemented as
a reduction in K+ conductance to mimic the net ex-
citability changes induced by 5-HT2A activation. In
each simulation, a focal perturbation—a square wave
of 0.5 Hz amplitude and 50 ms duration—was de-
livered to the excitatory population of the caudal
right middle frontal gyrus, and PCI was computed
across 15 different noise seeds. As shown in Fig-
ure 4 (left panel), the simulated psychedelic condition
did not produce a statistically significant change in
PCI compared to baseline, consistent with empirical

TMS-EEG data in humans under psilocybin [32]. The
distribution of PCI values was highly overlapping
between the two conditions, indicating preserved
large-scale causal integration despite the simulated
neurochemical modulation. The robustness analysis
(Figure 4, right panel) further confirmed that this lack
of significant PCI change was stable across a wide
range of excitatory and inhibitory reversal potential
values, with most parameter combinations showing
no significant difference (“ns”) between conditions.
To ensure that the choice of stimulation amplitude
did not bias PCI estimates, we systematically com-
pared baseline and reduced gK conditions across a
range of stimulation amplitudes (0.25–1 Hz) and mul-
tiple seeds. These control analyses revealed no signif-
icant differences in PCI outcomes across stimulation
values, thereby confirming the robustness of our re-
sults to stimulation amplitude within this range (see
Supplementary Fig. S6). On this basis, we selected an
amplitude of 0.5 Hz for all subsequent simulations,
balancing consistency with prior work and numerical
stability.
Simulated 5-HT2A receptor activation accounts for the high complexity of brain activity during psychedelic states

<figure><img src="image_6.png" /></figure>

**Figure 5: Spectral Power Analysis and Robustness to Parameter Variations.** Left) The plot shows seed- and channel-average spontaneous Power Spectral Density (PSD) plotted on a semi-log scale, highlighting the difference between the Baseline (black line) and reduced K+ conductance (blue line) conditions. The PSD was computed for both conditions and segmented into conventional EEG bands for analysis ($\delta, \theta, \alpha, \beta, \gamma$). The plot reveals a significant decrease in spectral power in the $\delta, \theta, \alpha$ bands under the reduced K+ leak conductance condition, indicated by the triple asterisks (t-test, *p* > 0.001). These results align with a known decrease in these bandwidths observed in literature under serotoninergic psychedelics. Shaded areas represent the standard deviation across simulations and channels, while solid lines represent the average PSD. Right) These heatmaps illustrate the normalized differences in spectral power for the $\delta, \theta, \alpha$ and $\beta$ bands, demonstrating the robustness of the observed effect. No significant difference was found in the gamma bands in any of the tested configurations. The x- and y-axes represent variations in the inhibitory and excitatory ($EL_i$ and $EL_e$) reversal potentials, respectively (baseline $EL_e$ = -63 and $EL_i$ = -65). The highlighted cells in green correspond to the parameters used in the figure on the left. The color gradient indicates the magnitude and direction of the normalized spectral power difference, with red signifying a positive difference and blue a negative difference (t-tests, **: $p < 0.001$, **: $p < 0.01$, *: $p < 0.05$, ns: $p > 0.05$). This figure demonstrates that the decrease in the $\delta, \theta, \alpha$ channels is a robust finding that holds true across a range of reversal potential values for both excitatory and inhibitory neural populations.

## 2.5 Spectral Power in Simulated Baseline and Psychedelic States

Alterations in spontaneous spectral power are among the most consistent electrophysiological signatures of serotoninergic psychedelic states. Across EEG and MEG studies these compounds reliably reduce power in low- and mid-frequency bands—particularly in the alpha range—while often increasing spontaneous signal diversity (Table 1). This pattern has been interpreted as reflecting a breakdown of oscillatory synchrony in canonical resting-state networks, such as the default mode network, in line with theories of increased neural entropy under psychedelics [52].

To assess whether our model reproduces these empirical spectral changes, we computed the channel-averaged power spectral density (PSD) for two simulated conditions: (i) baseline wakefulness and (ii) a psychedelic state induced by reducing K+ leak conductance to mimic the excitability effects of 5-HT<sub>2A</sub> receptor activation. The PSD was estimated for each condition, averaged across 15 independent noise re-

alizations, and decomposed into conventional EEG frequency bands ($\delta, \theta, \alpha, \beta, \gamma$).

As shown in Figure 5 (left panel), the psychedelic simulation exhibited a significant ($p < 0.001$) reduction in spectral power within the $\delta, \theta$, and $\alpha$ bands, closely mirroring human EEG and MEG findings under psilocybin [34, 35, 32], Lysergic acid diethylamide (LSD) [51], and dimethyltryptamine (DMT) [30]. No significant changes were observed in the $\beta$ or $\gamma$ bands. Shaded areas represent variability across simulations, with solid lines denoting the mean PSD.

The robustness analysis (Figure 5, right panel) confirmed that the power reduction in $\delta, \theta$, and $\alpha$ bands was not dependent on a narrow set of model parameters: across a broad range of excitatory and inhibitory reversal potential values, the effect persisted, with most parameter combinations yielding statistically significant decreases (absence of “ns” labels).

These results show that our model captures an important feature of psychedelic neurophysiology: a robust reduction in low-frequency oscillatory power,
especially in the alpha band, co-occurring with increased spontaneous complexity (Section 2.3) and preserved perturbational complexity (Section 2.4). This convergence between simulated and empirical data supports the model's validity in linking receptor-level modulation to large-scale electrophysiological changes observed in human brain recordings.

## 3 Discussion

Understanding the intricate mechanisms by which serotonergic psychedelics alter human brain activity at the whole-brain level remains a significant challenge in neuroscience. The multiscale modeling framework presented here addresses this challenge from a modeling perspective by integrating biophysical details, from single-neuron properties to large-scale network dynamics, providing a principled approach to investigate the effects of these compounds. By explicitly incorporating the action of psychedelics on 5-HT₂A receptors and their influence on leak membrane K⁺ conductances, our model offers a mechanistic bridge between molecular-level perturbations and emergent macroscopic brain states.

A key finding of this work is the model's ability to robustly replicate the experimentally observed increase in brain complexity, as quantified by the Lempel-Ziv complexity (LZc) index, under simulated psychedelic conditions. This increase in LZc, characterized by asynchronous-irregular dynamics and enhanced complexity compared to the resting waking state, aligns with the "entropic brain hypothesis" [52]. Our results provide an important computational support for the notion that heightened neural signal diversity and information richness underpin the altered states of consciousness induced by these compounds.

Furthermore, the model accurately replicates experimental findings regarding the Perturbational Complexity Index (PCI), demonstrating that PCI increases, but not significantly under psychedelic drug administration. This crucial dissociation, where spontaneous complexity (LZc) is increased while perturbational complexity (PCI) is preserved, highlights the distinct nature of these two measures of consciousness. The model posits that psychedelics enhance local excitability without compromising the brain's integrative causal architecture under perturbation. This suggests that while the content and experience of consciousness are profoundly modulated, the underlying capacity for consciousness, as indexed by PCI, is maintained. Such findings underscore the domain-specific sensitivity of different complexity metrics and emphasize the necessity of employing complementary measures to fully characterize altered states of consciousness, including high-order information theoretic measures [65, 66] and dynamic and emergence frameworks [67, 68, 69, 70].

Our simulations further demonstrate that the marked reduction in low-frequency spectral power ($\delta, \theta, \alpha$) observed under diminished K⁺ leak conductance is a robust phenomenon, persisting across a wide range of excitatory and inhibitory reversal potentials. This effect resonates with empirical reports of reduced alpha and low-frequency power under serotonergic psychedelics, suggesting that alterations in intrinsic conductances may provide a mechanistic substrate for these characteristic spectral signatures.

Beyond the specific context of psychedelics, this multiscale modeling framework represents a significant advancement towards the creation of "digital twins" of the brain [71]. By allowing for the systematic investigation of how changes at the molecular and cellular levels translate into macroscopic patterns of brain activity, the framework enhances the precision and applicability of whole-brain modeling for both fundamental neuroscience and potential clinical applications. Future work will involve further validation against a broader range of experimental data, exploring the effects of different receptor subtypes and their distributions, and investigating the framework's utility in modeling other pharmacological interventions and pathological brain states.

## 4 Current limitations and future research

While the present modeling framework provides a bottom-up bridge between receptor-level perturbations and emergent whole-brain dynamics, several limitations must be acknowledged. First, although our approach incorporates biologically realistic mechanisms such as 5-HT₂A receptor-mediated modulation of potassium leak conductance, the model currently focuses on a single receptor system. Psychedelic compounds act on multiple neurotransmitter receptors, including dopamine and glutamate [72, 73, 18], whose contributions remain to be explored in future extensions of this work. Incorporating these additional receptor targets, as well as their heterogeneous spatial distributions, would help to refine the generality and predictive power of the framework [74].

Second, our modeling approach relies on group-level structural connectivity data, which may not represent well individual differences in anatomy, receptor density, and pharmacokinetics. Integrating subject-specific connectomes and PET-based receptor maps could enable the development of personalized
Simulated 5-HT2A receptor activation accounts for the high complexity of brain activity during psychedelic states

<table><caption>**Table 1:** Summary of neurophysiological findings in psychedelic research. Studies are ordered chronologically. Arrows indicate: ↑ increase, ↓ decrease, ↔ no change, — not measured. LZ = Lempel-Ziv.</caption><thead><tr><th>Study</th><th>Method</th><th>Substance(s)</th><th>PCI</th><th>LZ Complexity</th><th>Spectral Power</th></tr></thead><tbody><tr><td>Muthukumaraswamy et al. (2013) [35]</td><td>MEG (spontaneous)</td><td>Psilocybin</td><td>—</td><td>↑</td><td>↓ Delta, Theta, Alpha (heterogeneous)</td></tr><tr><td>Carhart-Harris et al. (2016) [51]</td><td>MEG (spontaneous)</td><td>LSD</td><td>—</td><td>↑</td><td>↓ Alpha (heterogeneous)</td></tr><tr><td>Schartner et al. (2017) [51]</td><td>MEG (spontaneous)</td><td>Ketamine, LSD, Psilocybin</td><td>—</td><td>↑</td><td>↓ Alpha</td></tr><tr><td>Timmermann et al. (2019) [30]</td><td>EEG (spontaneous)</td><td>DMT</td><td>—</td><td>↑</td><td>↓ Alpha, Beta</td></tr><tr><td>Farnes et al. (2020) [31]</td><td>EEG (spontaneous) + TMS-EEG</td><td>Ketamine</td><td>↔</td><td>↑</td><td>↓ Theta, Alpha, and low Gamma (Eyes Closed)</td></tr><tr><td>Ort et al. (2023) [32]</td><td>EEG (spontaneous) + TMS-EEG</td><td>Psilocybin</td><td>↔</td><td>↑ (Eyes Closed)</td><td>↓ Theta, Alpha</td></tr><tr><td>Murray et al. (2024) [33]</td><td>EEG (spontaneous)</td><td>LSD</td><td>—</td><td>↑</td><td>↓ Delta, Theta</td></tr></tbody></table>

models of psychedelic action, providing new opportunities for linking emergent neurobiological features with individual variability. Such "digital twins" modeling approaches are particularly relevant for clinical translation.

Third, the current analysis emphasizes univariate measures of complexity and spectral properties. While these measures capture robust aspects of brain states induced by psychedelics, higher-order and dynamical functional interactions are likely essential to explain subjective phenomena such as ego dissolution and complex visual imagery. Incorporating recent developments in dynamic network analysis will be necessary to extend the explanatory scope of whole-brain models [75, 76, 67].

Future research should also explore non-linear dose-response relationships, as emerging evidence suggests that psychedelic effects are not strictly proportional to receptor occupancy and that threshold-like dynamics may govern the emergence of novel brain states. Extending the present framework to systematically investigate such non-linearities could reveal critical transitions in brain dynamics with direct relevance for dosing strategies.

# 5 Methods

We describe here the details of the multiscale whole-brain model: we start at the cellular level with the introduction of the single-cell model, then we move to the mesoscale with the construction of a neuronal network and the formulation of a mean-field model describing the mesoscopic network activity, and we end at the whole brain scale using a large network

of mean-fields interconnected using the anatomical connectome, step realized through the integration in The Virtual Brain simulator (Figure 2).

## 5.1 Single cell model

The dynamics of the single neuron are based on the adaptive exponential integrate-and-fire model (AdEx) and are described by the following equations [36]:

$$c_m \frac{dv}{dt} = g_L(E_L - v) + g_L \Delta e^{\frac{v-v_{thr}}{\Delta}} - w + I_{syn} \quad (7)$$

$$\frac{dw}{dt} = -\frac{1}{\tau_w}(a(v - E_L) - w) + b \sum_{t_{sp}} \delta(t - t_{sp}) \quad (8)$$

where $v$ is the membrane potential, $c_m = 200$ pF is the membrane capacitance and $g_L = 10$ nS (in the baseline condition) is the leak conductance. The leakage reversal potential $E_L$ equals -65 mV and -63 mV for inhibitory and excitatory cells respectively in the baseline. The exponential term $\Delta$ has a different strength for the two types of cells (2 mV for excitatory and 0.5 mV for inhibitory cells). $W$ and $I_{syn}$ describe the adaptation and synaptic current respectively. Parameter $a$ describes the sub-threshold adaptation and $b$ describes the spike-triggered adaptation. A spike is generated when the membrane potential exceeds a voltage threshold $v_{thr} = -50$ mV at time $t_{sp}$. The neuron's membrane potential is subsequently reset to a resting voltage $v_{rest} = -65$ mV and fixed to that value for a refractory period $T_{refr} = 5$ ms. The Dirac $\delta$-function indicates that whenever a neuron fires at time $t_{sp}(k)$, the adaptation current $W$ is incremented by an amount $b$. Based on physiological
characteristics [77], the inhibitory cells are modelled as fast-spiking (FS) neurons, exhibiting no adaptation (aᵢ = bᵢ = 0). Excitatory neurons are modelled as regular spiking (RS) cells, characterised by lower excitability due to the presence of spiking frequency adaptation (bₑ varies in our simulations, aₑ = 0nS, and the adaptation time constant τₑ = 500ms).

## 5.2 Synaptic model

Each neuron k receives a synaptic current I<sub>syn</sub>, which corresponds to the spiking activity of all presynaptic neurons j ∈ pre(k) of neuron k. I<sub>syn</sub> can be decomposed to the input received from excitatory (E) and inhibitory (I) presynaptic spikes, so as

$$I_{\text{syn}} = G_{\text{syn}}^{e}(E_e - v) + G_{\text{syn}}^{i}(E_i - v), \quad (9)$$

where E<sub>e</sub> = 0mV (E<sub>i</sub> = -80mV) is the excitatory (inhibitory) reversal potential and G<sub>syn</sub><sup>e</sup> (G<sub>syn</sub><sup>i</sup>) the excitatory (inhibitory) synaptic conductance. Synaptic conductances were modeled by a decaying exponential function that sharply increases by a fixed amount Q<sub>x</sub>, at each spiking time (t<sub>sp,j<sub>x</sub></sub>) of a presynaptic neuron j<sub>x</sub>:

$$G_{\text{syn}}^{x}(t) = Q_{x} \sum_{t_{\text{sp,jx}}} \Theta(t - t_{\text{sp,jx}}) e^{-\frac{t - t_{\text{sp,jx}}}{\tau_{x}}}, \quad (10)$$

where x corresponds to the population type (x ∈ {e, i}), Θ is the Heaviside function, τ<sub>x</sub> is the characteristic decay time of synaptic conductances (varied in our simulations), and Q<sub>x</sub> is the quantal conductance (Q<sub>e</sub> = 1.5nS, Q<sub>i</sub> = 5nS). The sum runs over all the spiking times of excitatory (t<sub>sp,j<sub>e</sub></sub>) or inhibitory (t<sub>sp,j<sub>i</sub></sub>) presynaptic neurons j<sub>e</sub> or j<sub>i</sub> of neuron k.

## 5.3 Spiking network model

In this work we considered a network representing a prototypical cortical circuit, specifically a single cortical column [42]. We studied a network of N = 10⁴ exponential-integrate-and-fire neurons that displayed spike-frequency adaptation [77, 44]. The neurons were connected over a topologically random network with a probability of connection between two neurons equal to p = 5%. The network was composed of two populations of inhibitory and excitatory neurons, with the inhibitory neurons consisting the 20% of the whole network size.

## 5.4 Mean-field model

The mean-field used here is based on a bottom-up formalism derived from a Master Equation [43], and which has recently been extended to include the effects of adaptation [44, 38]. The mean-field equations of the system can be written as:

$$T \frac{\partial v_{\mu}}{\partial t} = (F_{\mu} - v_{\mu}) + \frac{1}{2} c_{\lambda\eta} \frac{\partial^2 F_{\mu}}{\partial v_{\lambda} \partial v_{\eta}} \quad (11)$$

$$T \frac{\partial c_{\lambda\eta}}{\partial t} = \delta_{\lambda\eta} \frac{F_{\lambda}(1/T - F_{\eta})}{N_{\lambda}} + (F_{\lambda} - v_{\lambda})(F_{\eta} - v_{\eta}) \\ + \frac{\partial F_{\lambda}}{\partial v_{\mu}} c_{\eta\mu} + \frac{\partial F_{\eta}}{\partial v_{\mu}} c_{\lambda\mu} - 2c_{\lambda\eta} \quad (12)$$

$$\frac{\partial W}{\partial t} = \frac{1}{\tau_w} (a(\mu_V(v_e, v_i, W) - E_L) - W) + bu(13)$$

where μ = {e, i} is the neural population index (excitatory or inhibitory), v<sub>μ</sub> the mean firing rate of the corresponding population, c<sub>λη</sub> the covariance between populations (λ and η), W is the mean adaptation, and T is the mean-field characteristic time constant.

The function F<sub>μ</sub> = F<sub>μ</sub>(v<sub>e</sub>, v<sub>i</sub>, W) is the transfer function (TF) of a neuron of type μ, i.e. its output firing rate when receiving inhibitory and excitatory inputs with rates v<sub>e</sub> and v<sub>i</sub>, and adaptation level W. The TF can be derived following a semi-analytic approach [44], where the output firing rate of a neuron can be written as a function of its mean subthreshold membrane voltage μ<sub>V</sub>, its standard deviation σ<sub>V</sub>, and its time correlation decay time τ<sub>V</sub>:

$$F = v_{\text{out}} = \frac{1}{2\tau_V} \cdot \text{Erfc} \left( \frac{V_{\text{thr}}^{\text{eff}} - \mu_V}{\sqrt{2}\sigma_V} \right). \quad (14)$$

where (μ<sub>V</sub>, σ<sub>V</sub>, τ<sub>V</sub>) are calculated as a function of the input firing rates (v<sub>E</sub>, v<sub>I</sub>) and the adaptation intensity W following the equations described in [38]. Specifically, the mean membrane potential is calculated as the stationary solution under static conductances driven by the average synaptic input generated by firing rates (v<sub>E</sub>, v<sub>I</sub>). This input determines the mean (μ<sub>G<sub>e</sub></sub>, μ<sub>G<sub>i</sub></sub>) and standard deviation (σ<sub>G<sub>e</sub></sub>, σ<sub>G<sub>i</sub></sub>) of conductances for excitatory and inhibitory processes. Assuming Poissonian spike statistics (arising from asynchronous irregular dynamics), these values are expressed as:

$$\mu_{G_s}(v_e, v_i) = v_s K_s \tau_s Q_s, \quad (15)$$

$$\sigma_{G_s}(v_e, v_i) = \sqrt{v_s K_s \tau_s^2 Q_s}. \quad (16)$$

where $K_s = pN_s$ and $s = \{e, i\}$

The mean input conductance (μ<sub>G</sub>) and the effective membrane time constant (τ<sub>m</sub><sup>eff</sup>) are given by:

$$\mu_G(v_e, v_i) = \mu_{G_e} + \mu_{G_i} + g_L, \quad (17)$$

$$\tau_m^{\text{eff}}(v_e, v_i) = \frac{C_m}{\mu_G}. \quad (18)$$
Simulated 5-HT2A receptor activation accounts for the high complexity of brain activity during psychedelic states

The average membrane potential for a given adaptation current $w$ is:

$$\mu_V(v_e, v_i, w) = \frac{\mu_{G_e} E_e + \mu_{G_i} E_i + g_L E_L - w}{\mu_G}. \quad (19)$$

Note that in this work the terms $g_L$ and $E_L$ are respectively expressed by (3) and (4).

The standard deviation ($\sigma_V$) and time constant ($\tau_V$) of voltage fluctuations are:

$$\sigma_V(v_e, v_i) = \sqrt{\sum_s K_s v_s \frac{(U_s \cdot \tau_s)^2}{2(\tau_m^{\text{eff}} + \tau_s)}}, \quad (20)$$

$$\tau_V(v_e, v_i) = \frac{\sum_s (K_s v_s (U_s \cdot \tau_s)^2)}{\sum_s (K_s v_s \frac{(U_s \cdot \tau_s)^2}{\tau_m^{\text{eff}} + \tau_s})}. \quad (21)$$

where $s = \{e, i\}$ and $U_s = Q_s \mu_G (E_s - \mu_V)$.

The $V_{\text{thr}}^{\text{eff}}$ in Eq. 14 is the phenomenological spike threshold voltage taken as a second-order polynomial:

$$V_{\text{eff}}^{\text{thr}}(\mu_V, \sigma_V, \tau_V^N) = P_0 + \sum_{x \in \{\mu_V, \sigma_V, \tau_V^N\}} P_x \left( \frac{x - x_0}{\delta x_0} \right) + \sum_{x, y \in \{\mu_V, \sigma_V, \tau_V^N\}^2} P_{xy} \left( \frac{x - x_0}{\delta x_0} \right) \left( \frac{y - y_0}{\delta y_0} \right), \quad (22)$$

where $\tau_V^N = \frac{\tau_V G_i}{C_m}$ a non-dimensional quality. The polynomial coefficients of $P$ are determined through a fitting of the TF template to the output firing rate of individual neuron simulations, varying both inhibitory and excitatory inputs. The values for the normalization of the fluctuation regime were set following previous work [44, 38]: $\mu_V^0 = -60$ mV, $\sigma_V^0 = 0.004$ mV, $(\tau_V^N)^0 = 0.5$, $\delta\mu_V^0 = 0.001$ mV, $\delta c_V^0 = 0.006$ mV and $\delta(\tau_V^N)^0 = 1$. The fit is performed for each of the cell types included in the network, so in our case for fast and regular spiking cells. As shown in [38] the mean-field predictions work even far from the fitting point of the TF, so in all the examples showcased in this work we used previously calculated $P$, fitted for similar parameterizations of the two types of cells (see Table 2).

Taking the first-order equations by disregarding the covariance term dynamics in Equation 11, the mean-field model for the excitatory and inhibitory populations representing a cortical volume can be written as:

$$T \frac{dv_e}{dt} = \mathcal{F}_e(v_e + v_{\text{aff}}, v_i, W) - v_e \quad (23)$$

$$T \frac{dv_i}{dt} = \mathcal{F}_i(v_e + v_{\text{aff}}, v_i, W) - v_i \quad (24)$$

$$\frac{dW}{dt} = -\frac{W}{\tau_w} + b v_e + \frac{a}{\tau_w} (\mu_V(v_e, v_i, W) - E_L) \quad (25)$$

where $v_{\text{aff}}$ denotes an afferent input to both type of populations which we write as:

$$v_{\text{aff}}(t) = v_{\text{drive}} + \sigma \xi(t), \quad (26)$$

where $\sigma = 3.5$ and $\xi(t)$ denotes an Ornstein-Uhlenbeck (OU) process, of the form:

$$d\xi(t) = -\xi(t) \frac{dt}{\tau_{\text{OU}}} + dW_t \quad (27)$$

where $\tau_{\text{OU}} = 5$ms is the timescale of the OU process and $dW_t$ is a Wiener process of amplitude one and zero average. The $v_{\text{drive}}$ received different values in our simulations depending on the microscopic parameter that was selectively changed.

## 5.5 Networks of mean-field models

To model whole-brain dynamics, a network of mean-fields is defined, where each mean-field describes the activity of a brain region. The interactions between the mean-fields are excitatory while the inhibitory connections are preserved solely on the regional level. Extending on the single mean-field equations (Eq. 23), now the population activity of each region is given by the following equations:

$$T \frac{dv_e(k)}{dt} = \mathcal{F}_e[v_e^{\text{input}}(k), v_i(k), W(k)] - v_e(k) \quad (28)$$

$$T \frac{dv_i(k)}{dt} = \mathcal{F}_i[v_e^{\text{input}}(k), v_i(k), W(k)] - v_i(k) \quad (29)$$

$$\tau_w \frac{dW(k)}{dt} = -W(k) + b \tau_w v_e(k) \quad (30) \\ +a(\mu_V(v_e(k), v_i(k), W(k)) - E_L)$$

where $v_e(k)$ ($v_i(k)$) describes the firing rates of the excitatory (inhibitory) population of the region ($k$), $W(k)$ the adaptation of the population, and $v_e^{\text{input}}(k)$ is the total excitatory synaptic input that the region receives from the rest of the nodes of the network, given by the following equation:

$$v_e^{\text{input}}(k) = v_e(k) + v_{\text{aff}}(k) + G \sum_j C_{jk} v_e(j, t - \|j - k\| / v_c) \quad (31)$$

where the sum runs over all nodes $j$ connected to node $k$, and $C_{jk}$ the connection strength between $j$ and $k$ (and is equal to 1 for $j = k$), and the global coupling factor $G$ rescales the connection strength while maintaining the ratio between different values. The term $\|j - k\|$ is the distance between the nodes $j$ and $k$, while $v_c$ is the speed with which the signal propagates along the axis so that the model accounts for the delay of axonal propagation. Here $v_{\text{aff}}$ is defined as in Eq. 26.
<table><caption>**Table 2:** (Rounded) Fitting parameters expressed in mV</caption><thead><tr><th>Cell Type</th><th>P<sub>0</sub></th><th>P<sub>μV</sub></th><th>P<sub>σV</sub></th><th>P<sub>τ<sub>N</sub><sup>V</sup></sub></th><th>P<sub>μ<sub>V</sub><sup>2</sup></sub></th><th>P<sub>σ<sub>V</sub><sup>2</sup></sub></th><th>P<sub>τ<sub>N</sub><sup>V<sup>2</sup></sup></sub></th><th>P<sub>μ<sub>V</sub>σ<sub>V</sub></sub></th><th>P<sub>μ<sub>V</sub>τ<sub>N</sub><sup>V</sup></sub></th><th>P<sub>σ<sub>V</sub>τ<sub>N</sub></sub></th></tr></thead><tbody><tr><td>RS-cell</td><td>-49.83</td><td>5.064</td><td>-23.47</td><td>2.295</td><td>-0.411</td><td>10.55</td><td>-36.59</td><td>7.437</td><td>1.265</td><td>-40.72</td></tr><tr><td>FS-cell</td><td>-51.49</td><td>4.004</td><td>-8.352</td><td>0.241</td><td>-0.507</td><td>1.435</td><td>-14.69</td><td>4.503</td><td>2.847</td><td>-15.36</td></tr></tbody></table>

<table><caption>**Table 3:** Model Parameters.</caption><thead><tr><th>Parameter name</th><th>Symbol</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td colspan="4"><strong>Cellular Properties</strong></td></tr><tr><td>Sodium reversal potential</td><td>E<sub>Na</sub></td><td>50</td><td>mV</td></tr><tr><td>Potassium reversal potential</td><td>E<sub>K</sub></td><td>-90</td><td>mV</td></tr><tr><td>Leak conductance (baseline)</td><td>g<sub>L</sub></td><td>10</td><td>nS</td></tr><tr><td>Membrane capacitance</td><td>C<sub>m</sub></td><td>200</td><td>pF</td></tr><tr><td>Resting voltage</td><td>v<sub>rest</sub></td><td>-65</td><td>mV</td></tr><tr><td>Action Potential threshold</td><td>v<sub>thr</sub></td><td>-50</td><td>mV</td></tr><tr><td>Refractory period</td><td>T<sub>refr</sub></td><td>5</td><td>ms</td></tr><tr><td>Adaptation time constant</td><td>τ<sub>w</sub></td><td>500</td><td>ms</td></tr><tr><td colspan="4"><strong>Excitatory Neuron</strong></td></tr><tr><td>Leak reversal potential (baseline)</td><td>E<sub>L</sub></td><td>-63</td><td>mV</td></tr><tr><td>Spike sharpness</td><td>Δ<sub>e</sub></td><td>2</td><td>mV</td></tr><tr><td>Adaptation current increment</td><td>b<sub>e</sub></td><td>5</td><td>pA</td></tr><tr><td>Adaptation conductance</td><td>a<sub>e</sub></td><td>0</td><td>nS</td></tr><tr><td colspan="4"><strong>Inhibitory Neuron</strong></td></tr><tr><td>Leak reversal potential (baseline)</td><td>E<sub>L</sub></td><td>-65</td><td>mV</td></tr><tr><td>Spike sharpness</td><td>Δ<sub>i</sub></td><td>0.5</td><td>mV</td></tr><tr><td>Adaptation current increment</td><td>b<sub>i</sub></td><td>0</td><td>pA</td></tr><tr><td>Adaptation conductance</td><td>a<sub>i</sub></td><td>0</td><td>nS</td></tr><tr><td colspan="4"><strong>Synaptic Properties</strong></td></tr><tr><td colspan="4"><strong>Excitatory Neuron</strong></td></tr><tr><td>Reversal potential</td><td>E<sub>e</sub></td><td>0</td><td>mV</td></tr><tr><td>Quantal conductance</td><td>Q<sub>e</sub></td><td>1.5</td><td>nS</td></tr><tr><td>Decay time of synaptic conductance</td><td>τ<sub>e</sub></td><td>5</td><td>ms</td></tr><tr><td colspan="4"><strong>Inhibitory Neuron</strong></td></tr><tr><td>Reversal potential</td><td>E<sub>i</sub></td><td>-80</td><td>mV</td></tr><tr><td>Quantal conductance</td><td>Q<sub>i</sub></td><td>1.5</td><td>nS</td></tr><tr><td>Decay time of synaptic conductance</td><td>τ<sub>i</sub></td><td>5</td><td>ms</td></tr><tr><td colspan="4"><strong>Network Properties</strong></td></tr><tr><td>Total network size</td><td>N</td><td>10000</td><td></td></tr><tr><td>Global Coupling Factor</td><td>G</td><td>0.3</td><td></td></tr><tr><td>Signal Propagation Speed</td><td>V<sub>c</sub></td><td>4.0</td><td>mm/ms</td></tr><tr><td>Connectivity probability</td><td>p</td><td>0.05</td><td></td></tr><tr><td>Fraction of inhibitory cells</td><td>g<sub>ei</sub></td><td>0.2</td><td></td></tr></tbody></table>
Simulated 5-HT2A receptor activation accounts for the high complexity of brain activity during psychedelic states

A cortical parcellation of 68 regions was used. The connection strengths and tract lengths between the nodes were provided by the Berlin empirical data processing pipeline, based on human tractography methods [78]. To verify robustness of the results across connectomes, the same simulations were performed using a 412 regions connectome obtained from the cortical subset of regions of a 463 regions whole-brain connectome [79, 80]. The parameters used with this connectome were the same as with the 68 regions connectome, except for the targets $EL_e$ and $EL_i$, as smaller differences with baseline using this connectome appeared to be equivalent to larger differences in the 68 regions connectome. Similar results were obtained and can be found in the supplementary material (see Supplementary Figs. S5, S7, S8 and S9). We performed the simulations using The Virtual Brain (TVB) platform.

## 5.6 Heterogeneity

In the whole-brain model described the only parameters differentiating one region from another are found in the connectome, that is to say the weights between each region and their relative distance. However, brain regions differ in more varied ways, for instance in receptor density, cellular composition, etc. It therefore seems necessary for whole-brain models to take into account this heterogeneity between regions. We thus use a map of 5-HT₂A receptor densities [74, 81] to vary the activation and the effect on $g_K$ accordingly. Though this hypothesis could be refined with experimental data, we assumed a linear relationship between the density of 5-HT₂A receptors and the effect on the leak conductance. The leak conductance $g_{K_n}$ of region n is computed using its 5-HT₂A receptor density $d_n$ with a linear interpolation between $(0, g_{K_{baseline}})$ and $(\max_j(d_j), g_{K_{reduced}})$, such that an hypothetical region that would have a 5-HT₂A receptor density of 0 would not be affected by the reduced leak conductance and the region that has the highest 5-HT₂A receptor density has the lowest leak conductance (and hence the highest leak reversal potential). More formally:

$$g_{K_n} = g_{K_{baseline}} \left(1 - \frac{d_n}{\max_j(d_j)}\right) + g_{K_{reduced}} \frac{d_n}{\max_j(d_j)} \quad (32)$$

## 5.7 LZc calculation

Lempel-Ziv is a lossless compression algorithm that constructs a dictionary of substrings (words) sufficient to represent a string of characters given as input. A measure of complexity of the input string can then be obtained by simply counting the number

of words that were generated in the procedure, and can be thought of as the amount of distinct patterns necessary to summarize the data. In our case, our input data are excitatory firing rates obtained from each region of the connectome in the last second of a 5 seconds simulation. This represents 10000 data points for each region. As was done in previous studies for EEG and MEG signals [82], we apply the Hilbert transform on each channel, take the absolute, and binarize the data by defining as a threshold for each channel the mean value of the channel, such that each datapoint is converted to 1 if its value is over the mean and 0 otherwise. From this binarized data, we compute two different LZc-based measures, single-channel and multidimensional. Single-channel LZc consists in averaging the channels' individual LZc. The results of this measure on our simulations can be found in the supplementary material (Supplementary Figs. S4 and S5), and show an increase under K+ conductance reduction which holds for the two tested connectomes and multiple target leak reversal potential values, as for multidimensional LZc, which corresponds to experimental results. Multidimensional LZc (referred to as "LZc" in the paper) consists in first concatenating all binarized data to obtain a string of which the first n elements (with n the number of channels) are the binarized datapoints of each channel at the first timestep, the n following elements are the binarized datapoints of each channel at the second timestep, etc. We then apply the LZc algorithm on this binary string. For both these measures we apply a normalization scheme so that once LZc of a string is computed, we divide it by the LZc of the same string after a random permutation.

For statistical robustness we run 15 simulations with the same set of parameters but different random seeds and compute LZc each time.

## 5.8 Evoked activity and PCI calculation

The Perturbational Complexity Index (PCI), as proposed by Casali et al. [83] is the normalized Lempel-Ziv complexity of the evoked spatiotemporal patterns of cortical activation. This measure was proposed capture the complexity of the signal propagation evoked by an external stimulation. Here, the evoked potentials were generated using the whole-brain model, where a square wave stimulus was applied on the firing rates of the excitatory population on a single node (corresponding to the right premotor cortex). 15 different trials, each with a different stimulus onset and noise realization, were performed where the analysis concerns an interval of 300 ms pre- and post-stimulus. The firing rates of the excitatory populations were first normalized using the mean and standard deviation of the pre-stimulus in-
Simulated 5-HT2A receptor activation accounts for the high complexity of brain activity during psychedelic states

terval. To assess statistical significance, a bootstrap procedure was employed by repeatedly shuffling the pre-stimulus signals to generate a null distribution of surrogate values. A significance threshold was then determined from the null distribution. Finally, the normalized post-stimulus firing rates were compared to this threshold, resulting in a binary map of significant vectors $s(t)$ for each trial. The Lempel-Ziv complexity $LZ(S)$ was calculated for each of these vectors [9]. Finally the PCI is expressed as the ratio of $LZ(S)/H(S)$ where $H(S)$ is the spatial source entropy:

$$H(S) = -p(S = 0) \log_2(p(S = 0)) \quad (33) \\ - p(S = 1) \log_2(p(S = 1))$$

## 5.9 Power spectral density analysis

To quantify the frequency content of whole-brain activity, we computed the power spectral density (PSD) of excitatory population signals of all brain regions on the last second of 5 seconds simulations across 15 seeds for each condition. Let $x_r(t)$ denote the time series of firing rates of the excitatory population in region $r \in \{1, ..., N\}$, sampled at frequency $f_s = M/T$, where M is the number of time points and T the total recording duration. The discrete Fourier transform (DFT) of $x_r(t)$ is given by

$$X_r(k) = \sum_{m=0}^{M-1} x_r(m) e^{-i2\pi km/M}, \quad k = 0, ..., M-1,$$

with corresponding frequencies

$$f_k = \frac{kf_s}{M}, \quad f_k \in [0, f_s).$$

The power spectrum of region $r$ is defined as the squared magnitude of the Fourier transform,

$$P_r(f_k) = |X_r(k)|^2.$$

We then averaged the spectra across regions to obtain the mean PSD:

$$\bar{P}(f_k) = \frac{1}{N} \sum_{r=1}^{N} P_r(f_k),$$

## Data availability

The connectome used in the simulations for the Desikan Killiany parcellation can be found at https://github.com/mariasacha/paper_pipeline hub/tree/master/TVB/tvb_model_reference/data/connectivity [9]. The receptor map used in this paper for the Desikan Killiany parcellation are available at

https://github.com/netneurolab/hansen_receptors [74]. The connectome and receptor maps for the Laussane 463 parcellation can be found in https://github.com/singlesp/energy-landscape/tree/v1.0.3/data [84].

## Code availability

The source code of all simulations shown in this article is available via GitHub at https://github.com/Computational-NeuroPSI/simulated-serotonergic-receptors-tvb. Our code makes use of The Virtual Brain library, available at https://www.thevirtualbrain.org/tvb/zwei (version 2.9).

## Acknowledgments

Research supported by the CNRS, the ANR (FLAG-ERA BrainAct project) and the European Union (Human Brain Project H2020-785907, H2020-945539, Virtual Brain Twin project 101137289, EBRAINS-2.0 project 101147319).

## Author Contributions

A.D. conceived the model, H.M. performed the analysis and simulations, R.C wrote the final manuscript. All authors contributed to the interpretation of the results. R.C and A.D supervised the project.

## Competing Interests

The authors declare no competing interests.

## References

[1] Kristoffer AA Andersen, Robin Carhart-Harris, David J Nutt, and David Erritzoe. Therapeutic effects of classic serotonergic psychedelics: A systematic review of modern-era clinical studies. *Acta Psychiatrica Scandinavica*, 143(2):101–118, 2021.

[2] Robin Carhart-Harris, Bruna Giribaldi, Rosalind Watts, Michelle Baker-Jones, Ashleigh Murphy-Beiner, Roberta Murphy, Jonny Martell, Allan Blemings, David Erritzoe, and David J Nutt. Trial of psilocybin versus escitalopram for depression. *New England Journal of Medicine*, 384(15):1402–1411, 2021.
[3] Pim B van der Meer, Juan J Fuentes, Ad A Kaptein, Jan W Schoones, Marleen M De Waal, Anneke E Goudriaan, Kees Kramers, Arnt Schellekens, Metten Somers, Matthijs G Bossong, et al. Therapeutic effect of psilocybin in addiction: A systematic review. *Frontiers in Psychiatry*, 14:1134454, 2023.

[4] Lucas Arbabyazd, Kelly Shen, Zheng Wang, Martin Hofmann-Apitius, Petra Ritter, Anthony R McIntosh, Demian Battaglia, and Viktor Jirsa. Virtual connectomic datasets in alzheimer's disease and aging using whole-brain network dynamics modelling. *Eneuro*, 8(4), 2021.

[5] Leon Stefanovski, Jil Mona Meier, Roopa Kalsank Pai, Paul Triebkorn, Tristram Lett, Leon Martin, Konstantin Bülau, Martin Hofmann-Apitius, Ana Solodkin, Anthony Randal McIntosh, et al. Bridging scales in alzheimer's disease: Biological framework for brain simulation with the virtual brain. *Frontiers in Neuroinformatics*, 15:630172, 2021.

[6] Viktor Jirsa, Huifang Wang, Paul Triebkorn, Meysam Hashemi, Jayant Jha, Jorge Gonzalez-Martinez, Maxime Guye, Julia Makhalova, and Fabrice Bartolomei. Personalised virtual brain models in epilepsy. *The Lancet Neurology*, 22(5):443-454, 2023.

[7] Andrea I Luppi, Joana Cabral, Rodrigo Cofre, Pedro AM Mediano, Fernando E Rosas, Abid Y Qureshi, Amy Kuceyeski, Enzo Tagliazucchi, Federico Raimondo, Gustavo Deco, et al. Computational modelling in disorders of consciousness: Closing the gap towards personalised models for restoring consciousness. *NeuroImagage*, 275:120162, 2023.

[8] Katrin Amunts, Markus Axer, Swati Banerjee, Lise Bitsch, Jan G Bjaalie, Philipp Brauner, Andrea Brovelli, Navona Calarco, Marcel Carrere, Svenja Caspers, et al. The coming decade of digital brain research: A vision for neuroscience at the intersection of technology and computing. *Imaging Neuroscience*, 2024.

[9] Maria Sacha, Federico Tesler, Rodrigo Cofre, and Alain Destexhe. A computational approach to evaluate how molecular mechanisms impact large-scale brain activity. *Nature Computational Science*, 5(5):405-417, 2025.

[10] Gustavo Deco, Josephine Cruzat, Joana Cabral, Gitte M Knudsen, Robin L Carhart-Harris, Peter C Whybrow, Nikos K Logothetis, and

Morten L Kringelbach. Whole-brain multimodal neuroimaging model using serotonin receptor maps explains non-linear functional effects of lsd. *Current biology*, 28(19):3065-3074, 2018.

[11] Beatrice M Jobst, Selen Atasoy, Adrián Ponce-Alvarez, Ana Sanjuán, Leor Roseman, Mendel Kaelen, Robin Carhart-Harris, Morten L Kringelbach, and Gustavo Deco. Increased sensitivity to strong perturbations in a whole-brain model of lsd. *Neuroimage*, 230:117809, 2021.

[12] Rubén Herzog, Pedro AM Mediano, Fernando E Rosas, Paul Lodder, Robin Carhart-Harris, Yonatan Sanz Perl, Enzo Tagliazucchi, and Rodrigo Cofre. A whole-brain model of the neural entropy increase elicited by psychedelic drugs. *Scientific reports*, 13(1):6244, 2023.

[13] I Mindlin, R Herzog, L Belloli, D Manasova, M Monge-Asensio, J Vohryzek, A Escrichs, Naji Alnagger, P Núñez, Olivia Gosseries, et al. Whole brain modelling for simulating pharmacological interventions on patients with disorders of consciousness. *Communications Biology*, 7(1):1176, 2024.

[14] Vincent Beliveau, Melanie Ganz, Ling Feng, Brice Ozenne, Liselotte Højgaard, Patrick M Fisher, Claus Svarer, Douglas N Greve, and Gitte M Knudsen. A high-resolution in vivo atlas of the human brain's serotonin system. *Journal of Neuroscience*, 37(1):120-128, 2017.

[15] Rodrigo Cofré, Rubén Herzog, Pedro AM Mediano, Juan Piccinini, Fernando E Rosas, Yonatan Sanz Perl, and Enzo Tagliazucchi. Whole-brain models to explore altered states of consciousness from the bottom up. *Brain Sciences*, 10(9):626, 2020.

[16] Clinton E Canal. Serotonergic psychedelic: experimental approaches for assessing mechanisms of action. *New Psychoactive Substances: Pharmacology, Clinical, Forensic and Analytical Toxicology*, pages 227-260, 2018.

[17] Samuel T Slocum, Jeffrey F DiBerto, and Bryan L Roth. Molecular insights into psychedelic drug action. *Journal of neurochemistry*, 162(1):24-38, 2022.

[18] Alex C Kwan, David E Olson, Katrin H Preller, and Bryan L Roth. The neural basis of psychedelic action. *Nature Neuroscience*, 25(11):1407-1419, 2022.

[19] Alaina M Jaster and Javier González-Maeso. Mechanisms and molecular targets surrounding
the potential therapeutic effects of psychedelics.
*Molecular psychiatry*, 28(9):3595-3612, 2023.

[20] Paula Sanz-Leon, Stuart A Knock, Andreas Spiegler, and Viktor K Jirsa. Mathematical framework for large-scale brain network modeling in the virtual brain. *Neuroimage*, 111:385-430, 2015.

[21] Jason Wallach, Andrew B. Cao, Maggie M. Calkins, Andrew J. Heim, Janelle K. Lanham, Emma M. Bonniwell, Joseph J. Hennessey, Hai-ley A. Bock, Emilie I. Anderson, Alexander M. Sherwood, Hamilton Morris, Robbin de Klein, Adam K. Klein, Bruna Cuccurazzu, James Gam-rat, Tilka Fannana, Randy Zauhar, Adam L. Hal-berstadt, and John D. McCorvy. Identification of 5-ht2a receptor signaling pathways associated with psychedelic potential. *Nature Communica-tions*, 14(1):8221, Dec 2023.

[22] G. K. Aghajanian and G. J. Marek. Serotonin and hallucinogens. *Neuropsychopharmacology*, 21(1):16-23, Aug 1999.

[23] C P VanderMaelen and G K Aghajanian. Intra-cellular studies showing modulation of facial motoneuron excitability by serotonin. *Nature*, 287(5780):346-347, September 1980.

[24] C P Vandermaelen and G K Aghajanian. Serotonin-induced depolarization of rat facial motoneurons in vivo: comparison with amino acid transmitters. *Brain Res*, 239(1):139-152, May 1982.

[25] R C Ma and N J Dun. Excitation of lateral horn neurons of the neonatal rat spinal cord by 5-hydroxytryptamine. *Brain Res*, 389(1-2):89-98, January 1986.

[26] R Araneda and R Andrade. 5-hydroxytryptamine2 and 5-hydroxytryptamine1a receptors mediate opposing responses on membrane excitability in rat association cortex. *Neuroscience*, 40(2):399-412, 1991.

[27] David A Mccormick. Neurotransmitter actions in the thalamus and cerebral cortex and their role in neuromodulation of thalamocortical activity. *Progress in Neurobiology*, 39:337, 1992.

[28] Eiichiro Tanaka and R Alan North. Actions of 5-hydroxytryptamine on neurons of the rat cingulate cortex. Technical report, 1993.

[29] Michael M Schartner, Robin L Carhart-Harris, Adam B Barrett, Anil K Seth, and Suresh D Muthukumaraswamy. Increased spontaneous meg signal diversity for psychoactive doses of

ketamine, lsd and psilocybin. *Scientific reports*, 7(1):46421, 2017.

[30] Christopher Timmermann, Leor Roseman, Michael Schartner, Raphael Milliere, Lea T Williams, David Erritzoe, Suresh Muthukumaraswamy, Poppy Ashton, Charlotte Martial, Helena Cassol, Steven Laureys, David Nutt, and Robin Carhart-Harris. Neural correlates of the dmt experience assessed with multivariate eeg. *Scientific Reports*, 9(1):16324, 2019.

[31] Nadine Farnes, Bjørn E Juel, André S Nilsen, Luis G Romundstad, and Johan F Storm. Increased signal diversity/complexity of spontaneous eeg, but not evoked eeg responses, in ketamine-induced psychedelic state in humans. *Plos one*, 15(11):e0242056, 2020.

[32] Andres Ort, John W Smallridge, Simone Sarasso, Silvia Casarotto, Robin Von Rotz, Andrea Casanova, Erich Seifritz, Katrin H Preller, Giulio Tononi, and Franz X Vollenweider. Tms-eeg and resting-state eeg applied to altered states of consciousness: oscillations, complexity, and phenomenology. *IsScience*, 26(5), 2023.

[33] Conor H Murray, Joel Frohlich, Connor J Hag-garty, Ilaria Tare, Royce Lee, and Harriet de Wit. Neural complexity is increased after low doses of lsd, but not moderate to high doses of oral thc or methamphetamine. *Neuropsychopharmacology*, 49(7):1120-1128, 2024.

[34] Michael Kometer, Thomas Pokorny, Erich Seifritz, and Franz X Volleinweider. Psilocybin-induced spiritual experiences and insightfulness are associated with synchronization of neuronal oscillations. *Psychopharmacology (Berl)*, 232(19):3663-3676, August 2015.

[35] Suresh D. Muthukumaraswamy, Robin L. Carhart-Harris, Rosalyn J. Moran, Matthew J. Brookes, Tim M. Williams, David Erritzoe, Ben Sessa, Andreas Papadopoulos, Mark Bolstridge, Krish D. Singh, Amanda Feilding, Karl J. Friston, and David J. Nutt. Broadband cortical desynchronization underlies the human psychedelic state. *Journal of Neuroscience*, 33(38):15171-15183, 2013.

[36] Romain Brette and Wolfram Gerstner. Adaptive exponential integrate-and-fire model as an effective description of neuronal activity. *Journal of neurophysiology*, 94(5):3637-3642, 2005.

[37] Richard Naud, Nicolas Marcille, Claudia Clopath, and Wolfram Gerstner. Firing patterns
Simulated 5-HT2A receptor activation accounts for the high complexity of brain activity during psychedelic states

in the adaptive exponential integrate-and-fire model. *Biological cybernetics*, 99:335-347, 2008.

[38] Matteo Di Volo, Alberto Romagnoni, Cristiano Capone, and Alain Destexhe. Biologically realistic mean-field models of conductance-based networks of spiking neurons with adaptation. *Neural computation*, 31(4):653-680, 2019.

[39] Jorin Overwiening, Federico Tesler, Domenico Guarino, and Alain Destexhe. A multiscale study of thalamic state-dependent responsiveness. *PLOS Computational Biology*, 20(12):e1012262, 2024.

[40] Federico Tesler, Alexander Kozlov, Sten Grillner, and Alain Destexhe. A multiscale model of striatum microcircuit dynamics. *bioRxiv*, pages 2023-12, 2023.

[41] Roberta Maria Lorenzi, Alice Geminiani, Yann Zerlaut, Marialaura De Grazia, Alain Destexhe, Claudia AM Gandini Wheeler-Kingshott, Fulvia Palesi, Claudia Casellato, and Egidio D'Angelo. A multi-layer mean-field model of the cerebellum embedding microstructure and population-specific dynamics. *PLOS Computational Biology*, 19(9):e1011434, 2023.

[42] Valentino Braitenberg and Almut Schüz. *Context: statistics and geometry of neuronal connectivity*. Springer Science & Business Media, 2013.

[43] Sami El Boustani and Alain Destexhe. A master equation formalism for macroscopic modeling of asynchronous irregular activity states. *Neural computation*, 21(1):46-100, 2009.

[44] Yann Zerlaut, Sandrine Chemla, Frederic Chavane, and Alain Destexhe. Modeling mesoscopic cortical dynamics using a mean-field model of conductance-based networks of adaptive exponential integrate-and-fire neurons. *Journal of Computational Neuroscience*, 44(1):45-61, 2018.

[45] Mallory Carlu, Omar Chehab, Leonardo Dalla Porta, Damien Depannemaecker, Charlotte Héricé, Maciej Jedynak, E Köksal Ersöz, Paolo Muratore, Selma Souihel, Cristiano Capone, et al. A mean-field approach to the dynamics of networks of complex neurons, from nonlinear integrate-and-fire to hodgkin-huxley models. *Journal of neurophysiology*, 123(3):1042-1051, 2020.

[46] Christoffer G Alexandersen, Chloé Duprat, Aitakin Ezzati, Pierre Houzelstein, Ambre

Ledoux, Yuhong Liu, Sandra Saghir, Alain Destexhe, Federico Tesler, and Damien Depannemaecker. A mean field to capture asynchronous irregular dynamics of conductance-based networks of adaptive quadratic integrate-and-fire neuron models. *Neural Computation*, pages 1-16, 2024.

[47] Rahul S Desikan, Florent Ségonne, Bruce Fischl, Brian T Quinn, Bradford C Dickerson, Debrah Blacker, Randy L Buckner, Anders M Dale, R Paul Maguire, Bradley T Hyman, et al. An automated labeling system for subdividing the human cerebral cortex on mri scans into gyral based regions of interest. *Neuroimage*, 31(3):968-980, 2006.

[48] Geraint Rees, Gabriel Kreiman, and Christof Koch. Neural correlates of consciousness in humans. *Nature Reviews Neuroscience*, 3(4):261-270, 2002.

[49] Christof Koch, Marcello Massimini, Melanie Boly, and Giulio Tononi. Neural correlates of consciousness: progress and problems. *Nature reviews neuroscience*, 17(5):307-321, 2016.

[50] Jacob Ziv and Abraham Lempel. Compression of individual sequences via variable-rate coding. *IEEE transactions on Information Theory*, 24(5):530-536, 2003.

[51] Michael M. Schartner, Andrea Pigorini, A. Gibbs, Gabriele Arnulfo, Simone Sarasso, Lionel Barrett, Lino Nobili, Marcello Massimini, Anil K. Seth, and Adam B. Barrett. Increased spontaneous meg signal diversity for psychoactive doses of ketamine, lsd, and psilocybin. *Scientific Reports*, 7:46421, 2017.

[52] Robin L Carhart-Harris, Robert Leech, Peter J Hellyer, Murray Shanahan, Amanda Feilding, Enzo Tagliazucchi, Dante R Chialvo, and David Nutt. The entropic brain: a theory of conscious states informed by neuroimaging research with psychedelic drugs. *Frontiers in Human Neuroscience*, 8:20, 2014.

[53] Robin L Carhart-Harris. The entropic brain: revisited. *Neuropharmacology*, 142:167-178, 2018.

[54] Anthony G Hudetz, Xiping Liu, Siveshigan Pil-lay, Melanie Boly, and Giulio Tononi. Propo-fol anesthesia reduces lempel-ziv complexity of spontaneous brain activity in rats. *Neuroscience letters*, 628:132-135, 2016.
Simulated 5-HT2A receptor activation accounts for the high complexity of brain activity during psychedelic states

[55] X-S Zhang, Rob J Roy, and Erik W Jensen. Eeg complexity as a measure of depth of anesthesia for patients. *IEEE transactions on biomedical engineering*, 48(12):1424–1433, 2001.

[56] Rain Ferenets, Ann Vanluchene, Tarmo Lipping, Björn Heyse, and Michel MRF Struys. Behavior of entropy/complexity measures of the electroencephalogram during propofol-induced sedation: dose-dependent effects of remifentanil. *Anesthesiology*, 106(4):696–706, 2007.

[57] Daniela Biggs, Gonzalo Boncompte, Juan C Pedemonte, Carlos Fuentes, and Luis I Cortínez. The effect of age on electroencephalogram measures of anesthesia hypnosis: A comparison of bis, alpha power, lempel-ziv complexity and permutation entropy during propofol induction. *Frontiers in Aging Neuroscience*, 14:910886, 2022.

[58] Gonzalo Boncompte, Vicente Medel, Luis I Cortínez, and Tomás Ossandón. Brain activity complexity has a nonlinear relation to the level of propofol sedation. *British Journal of Anaesthesia*, 127(2):254–263, 2021.

[59] Thomas Andrillon, Andreas Trier Poulsen, Lars Kai Hansen, Damien Léger, and Sid Kouider. Neural markers of responsiveness to the environment in human sleep. *Journal of Neuroscience*, 36(24):6583–6596, 2016.

[60] Arnfinn Aamodt, André Sevenius Nilsen, Rune Markhus, Anikó Kusztor, Fatemeh Hasanzadeh Moghadam, Nils Kauppi, Benjamin Thürer, Johan Frederik Storm, and Bjørn Erik Juel. Eeg lempel-ziv complexity varies with sleep stage, but does not seem to track dream experience. *Frontiers in Human Neuroscience*, 16:987714, 2023.

[61] Michael M Schartner, Andrea Pigorini, Steve A Gibbs, Gabriele Arnulfo, Simone Sarasso, Lional Barnett, Lino Nobili, Marcello Massimini, Anil K Seth, and Adam B Barrett. Global and local complexity of intracranial eeg decreases during nrem sleep. *Neuroscience of consciousness*, 2017(1):niw022, 2017.

[62] Daniel Abásolo, Samantha Simons, Rita Morgado da Silva, Giulio Tononi, and Vladyslav V Vyazovskiy. Lempel-ziv complexity of cortical activity during sleep and waking in rats. *Journal of neurophysiology*, 113(7):2742–2752, 2015.

[63] Adrien G. Casali, Olivia Gosseries, Massimini Rosanova, et al. A theoretically based index of consciousness independent of sensory processing and behavior. *Science Translational Medicine*, 5(198):198ra105, 2013.

[64] Stefano Casarotto, Angela Comanducci, Massimini Rosanova, et al. Stratification of unresponsive patients by an independently validated index of brain complexity. *Annals of Neurology*, 80(6):718–729, 2016.

[65] Andrea I Luppi, Pedro AM Mediano, Fernando E Rosas, Negin Holland, Tim D Fryer, John T O'Brien, James B Rowe, David K Menon, Daniel Bor, and Emmanuel A Stamatakis. A synergistic core for human brain evolution and cognition. *Nature neuroscience*, 25(6):771–782, 2022.

[66] Laouen Belloli, Pedro Mediano, Rodrigo Cofré, Diego Fernandez Slezak, and Rubén Herzog. Thoi: An efficient and accessible library for computing higher-order interactions enhanced by batch-processing. *arXiv preprint arXiv:2501.03381*, 2025.

[67] Pablo Castro, Andrea Luppi, Enzo Tagliazucchi, Yonatan S Perl, Lorina Naci, Adrian M Owen, Jacobo D Sitt, Alain Destexhe, and Rodrigo Cofré. Dynamical structure-function correlations provide robust and generalizable signatures of consciousness in humans. *Communications Biology*, 7(1):1224, 2024.

[68] Andrea I Luppi, Lynn Uhrig, Jordy Tasserie, Camilo M Signorelli, Emmanuel A Stamatakis, Alain Destexhe, Bechir Jarraya, and Rodrigo Cofre. Local orchestration of distributed functional patterns supporting loss and restoration of consciousness in the primate brain. *Nature communications*, 15(1):2171, 2024.

[69] Borjan Milinkovic, Lionel Barnett, Olivia Carter, Anil K Seth, and Thomas Andrillon. Capturing the emergent dynamical structure in biophysical neural models. *PLOS Computational Biology*, 21(5):e1012572, 2025.

[70] Fernando E Rosas, Pedro AM Mediano, Henrik J Jensen, Anil K Seth, Adam B Barrett, Robin L Carhart-Harris, and Daniel Bor. Reconciling emergences: An information-theoretic approach to identify causal emergence in multivariate data. *PLoS computational biology*, 16(12):e1008289, 2020.

[71] Huifang E Wang, Paul Triebkorn, Martin Breyton, Borana Dollomaja, Jean-Didier Lemarechal, Spase Petkoski, Pierpaolo Sorrentino, Damien Depannemaecker, Meysam Hashemi, and Viktor K Jirsa. Virtual brain twins: from basic neuroscience to clinical use. *National Science Review*, 11(5):nwae079, 2024.
Simulated 5-HT2A receptor activation accounts for the high complexity of brain activity during psychedelic states

[72] Pamela A Pierce and Stephen J Peroutka. Hal-lucinogenic drug interactions with neurotrans-mitter receptor binding sites in human cortex. *Psychopharmacology*, 97(1):118-122, 1989.

[73] NL Mason, KPC Kuypers, F Müller, J Reck-weg, DHY Tse, SW Toennes, NRPW Hutten, JFA Jansen, P Stiers, A Feilding, et al. Me, my-self, bye: regional alterations in glutamate and the experience of ego dissolution with psilocy-bin. *Neuropsychopharmacology*, 45(12):2003-2011, 2020.

[74] Justine Y Hansen, Golia Shafiei, Ross D Markello, Kelly Smart, Sylvia ML Cox, Mar-tin Nørgaard, Vincent Beliveau, Yanjun Wu, Jean-Dominique Gallezot, Étienne Aumont, et al. Mapping neurotransmitter systems to the struc-tural and functional organization of the human neocortex. *Nature neuroscience*, 25(11):1569-1581, 2022.

[75] Sepehr Mortaheb, Larry D Fort, Natasha L Mas-on, Pablo Mallaroni, Johannes G Ramaekers, and Athena Demertzi. Dynamic functional hy-perconnectivity after psilocybin intake is pri-marily associated with oceanic boundlessness. *Biological Psychiatry: Cognitive Neuroscience and Neuroimaging*, 9(7):681-692, 2024.

[76] Athena Demertzi, Enzo Tagliazucchi, Stanis-las Dehaene, Gustavo Deco, Pablo Barttfeld, Federico Raimondo, Charlotte Martial, Davinia Fernández-Espejo, Benjamin Rohaut, HU Voss, et al. Human consciousness is supported by dynamic complex patterns of brain signal coor-dination. *Science advances*, 5(2):eaat7603, 2019.

[77] Alain Destexhe. Self-sustained asynchronous irregular states and up-down states in thala-mic, cortical and thalamocortical networks of nonlinear integrate-and-fire neurons. *Journal of Computational Neuroscience*, 27(3):493, 2009.

[78] Michael Schirner, Simon Rothmeier, Viktor K Jirsa, Anthony Randal McIntosh, and Petra Ritter. An automated pipeline for constructing personalized virtual brains from multimodal neuroimaging data. *NeuroImage*, 117:343-357, 2015.

[79] Leila Cammoun, Xavier Gigandet, Djalel Meskaldji, Jean Philippe Thiran, Olaf Sporns, Kim Q Do, Philippe Maeder, Reto Meuli, and Patric Hagmann. Mapping the human connec-tome at multiple scales with diffusion spectrum mri. *Journal of neuroscience methods*, 203(2):386-397, 2012.

[80] Andrea I Luppi and Emmanuel A Stamatakis. Combining network topology and information theory to construct representative brain networks. *Network Neuroscience*, 5(1):96-124, 2021.

[81] Vincent Beliveau, Melanie Ganz, Ling Feng, Brice Ozenne, Liselotte Højgaard, Patrick M Fisher, Claus Svarer, Douglas N Greve, and Gitte M Knudsen. A high-resolution *in vivo* atlas of the human brain's serotonin system. *J. Neurosci.*, 37(1):120-128, January 2017.

[82] Michael Schartner, Robin Carhart-Harris, Adam B. Barrett, Anil K. Seth, and Suresh Muthukumaraswamy. Complexity of multi-dimensional spontaneous eeg decreases during propofol induced general anaesthesia. *PLoS ONE*, 10(8):e0133532, 2015.

[83] Adenauer G Casali, Olivia Gosseries, Mario Rosanova, Mélanie Boly, Simone Sarasso, Karina R Casali, Silvia Casarotto, Marie-Aurélie Bruno, Steven Laureys, Giulio Tononi, et al. A theoretically based index of consciousness independent of sensory processing and behavior. *Science translational medicine*, 5(198):198ra105-198ra105, 2013.

[84] S Parker Singleton, Andrea I Luppi, Robin L Carhart-Harris, Josephine Cruzat, Leor Roseman, David J Nutt, Gustavo Deco, Morten L Kringelbach, Emmanuel A Stamatakis, and Amy Kuceyeski. Receptor-informed network control theory links lsd and psilocybin to a flattening of the brain's control energy landscape. *Nature communications*, 13(1):5812, 2022.
Supplementary material for: Simulated serotonergic receptor stimulation reproduces the high complexity of psychedelic states at the whole-brain level

<figure><img src="image_7.png" /><figcaption>Figure S1: Increased firing rates with reduced K+ leak conductance. Average firing rates of excitatory populations are shown for the baseline condition (dark gray) and for the condition with reduced K+ leak conductance (light blue) where $EL_e = -61.2$ and $EL_i = -64.4$. For illustration, results are displayed for the 5th simulated second of spontaneous activity (seed 0), but the effect is robust across simulations and independent of the seed. Thin traces represent individual population firing rates, while bold lines indicate their mean dynamics over time.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure S2: 5-HT<sub>2A</sub> Density of receptors. The image displays a positron emission tomography (PET) derived map of serotonin 5-HT<sub>2A</sub> receptor density on the cortical surface using the Desikan-Killiany with 68 cortical brain areas. The color bar indicates the receptor density, with a continuous spectrum from minimum (Min) to maximum (Max) values. High receptor density is observed in regions such as the frontal cortex and posterior cingulate cortex.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure S3: **Effect of receptor maps on model dynamics.** The boxplots show the distribution of the LZc across regions for different receptor maps. The baseline condition corresponds to no reduction in *g*<sub>*K*</sub>. The 5-HT<sub>2A</sub> receptor condition implements a region-specific reduction in potassium conductance proportional to the empirical 5-HT<sub>2A</sub> receptor density map (as in the main simulations). Additional conditions tested the specificity of this effect: other receptor maps, a spatially randomized (shuffled) 5-HT<sub>2A</sub> map, uniformly distributed receptor densities sampled between the minimum and maximum values of the empirical 5-HT<sub>2A</sub> map (uniform), and a homogeneous map in which all regions were assigned the maximum receptor density (homogeneous). The asterisks indicate the level of statistical significance (t-tests, **: *p* < 0.001, **: *p* < 0.01, *: *p* < 0.05). These results show that the specific distribution of the 5-HT<sub>2A</sub> receptor map is not necessary to obtain a significative increase in LZc in our model.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure S4: Lempel-Ziv complexity for single channels (LZc single) and Robustness Analysis. Left) Difference in LZc single between the Baseline and K+ conductance reduction conditions. The box plot displays the median (red line), interquartile range (box), and whiskers extending to 1.5 times the interquartile range. The violin plot shows the data distribution, and individual data points are overlaid. A significant difference, as indicated by the three asterisks (t-test, *p* < 0.001), is observed between the two conditions, with reduced K+ leak conductance condition showing a higher LZc single value. Right) The heatmap demonstrates the robustness of the LZc results to changes in the target excitatory and inhibitory reversal potentials (baseline *EL*<sub>e</sub> = -63 and *EL*<sub>i</sub> = -65). The highlighted cell in green corresponds to the parameters used in the figure on the left. The color bar on the right shows the normalized LZc single difference, where red indicates a positive difference (increase in LZc single) and blue a negative difference. The asterisks indicate the level of statistical significance (t-tests, **: *p* < 0.001, **: *p* < 0.01, *: *p* < 0.05, ns: *p* > 0.05). The prevalence of significant differences across a wide range of parameter values confirms that the observed increase in LZc single under the reduced K+ leak conductance condition is a robust finding.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure S5: Lempel-Ziv complexity for single channels (LZc single) and Robustness Analysis using the cortical-reduced Lausanne463 connectome. Left) This figure shows the difference in LZc single between the Baseline and K+ conductance reduction conditions. The box plot displays the median (red line), interquartile range (box), and whiskers extending to 1.5 times the interquartile range. The violin plot shows the data distribution, and individual data points are overlaid. A significant difference, as indicated by the three asterisks (t-test, *p* < 0.001), is observed between the two conditions, with the K+ conductance reduction condition showing a higher LZc single value. This is consistent with an increase in complexity reported in the literature under serotoninergic psychodelics. Right) The heatmap demonstrates the robustness of the LZc results to changes in the target excitatory and inhibitory reversal potentials (baseline *EL*<sub>e</sub> = -63 and *EL*<sub>i</sub> = -65). The highlighted cell corresponds to the parameters used in the figure on the left. The color bar on the right shows the normalized LZc single difference, where red indicates a positive difference (increase in LZc single) and blue a negative difference. The asterisks indicate the level of statistical significance (t-tests, **: *p* < 0.001, **: *p* < 0.01, *: *p* < 0.05, ns: *p* > 0.05). The prevalence of significant differences across a wide range of parameter values confirms that the observed increase in LZc single under the K+ conductance reduction condition is a robust finding.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure S6: **Effect of stimulation amplitude on Perturbational Complexity Index (PCI) Violin/box plots show the distribution of differences in PCI between baseline simulations (no reduction in potassium leak conductance, *g*<sub>K</sub>) and simulations with reduced *g*<sub>K</sub>, across multiple seeds, for a range of stimulation amplitudes (0.25–1 Hz). For each seed, the PCI difference was computed as the value under baseline minus the corresponding value under leak-reduced conditions at the given stimulation amplitude. Statistical comparisons indicate no significant differences (ns) across stimulation values, demonstrating that the choice of stimulation amplitude does not systematically influence PCI outcomes. Based on these results, a stimulation amplitude of 0.5 Hz was chosen for all subsequent analyses, ensuring comparability with previous work while avoiding instabilities (explosions) observed at higher excitatory drive (large *EL*<sub>e</sub> relative to *EL*<sub>i</sub>). This control confirms that PCI findings are robust to stimulation amplitude within the tested range.**</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure S7: Lempel-Ziv complexity (LZc) and Robustness Analysis using the cortical-reduced Lausanne463 connectome. Left) This figure shows the difference in Lempel-Ziv complexity (LZc) between the Baseline and K+ conductance reduction conditions. The box plot displays the median (red line), interquartile range (box), and whiskers extending to 1.5 times the interquartile range. The violin plot shows the data distribution, and individual data points are overlaid. A significant difference, as indicated by the three asterisks (t-test, *p* < 0.001), is observed between the two conditions, with the K+ conductance reduction condition showing a higher LZc value. This is consistent with an increase in complexity reported in the literature under serotoninergic psychedelics. Right) The heatmap demonstrates the robustness of the LZc results to changes in the target excitatory and inhibitory reversal potentials (baseline *EL*<sub>e</sub> = -63 and *EL*<sub>i</sub> = -65). The highlighted cell corresponds to the parameters used in the figure on the left. The color bar on the right shows the normalized LZc difference, where red indicates a positive difference (increase in LZc) and blue a negative difference. The asterisks indicate the level of statistical significance (t-tests, ***: *p* < 0.001, **: *p* < 0.01, *: *p* < 0.05, ns: *p* > 0.05). The prevalence of significant differences across a wide range of parameter values confirms that the observed increase in LZc under the K+ conductance reduction condition is a robust finding.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure S8: Perturbational Complexity Index (PCI) Analysis and Robustness using the cortical-reduced Lausanne463 connectome. Left) The box plot compares the Perturbational Complexity Index (PCI) values between the Baseline and K+ conductance reduction conditions. The central red lines represent the median PCI values, while the boxes indicate the interquartile range (IQR). The whiskers extend to 1.5 times the IQR. The violin plots show the distribution of the simulations and individual points are overlaid. The results show no statistically significant difference ("ns", t-test, *p* > 0.05) in PCI between the two conditions. Right) The heatmap demonstrates the robustness of the PCI results under changes in key model parameters. The x- and y-axes represent variations in the target excitatory and inhibitory reversal potentials. The highlighted cell corresponds to the parameters used in the figure on the left. The color scale represents the normalized PCI difference, with red indicating a positive difference and blue a negative difference. The abundance of "ns" (not significant) labels across the heatmap confirms that the lack of a significant change in PCI is a robust finding that remain true across a wide range of these parameter values.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure S9: Spectral Power Analysis and Robustness to Parameter Variations using the cortical-reduced Lausanne463 connectome. Left) The plot shows channel-average spontaneous Power Spectral Density (PSD) plotted on a semi-log scale, highlighting the difference between the Baseline (black line) and K+ conductance reduction (blue line) conditions. The PSD was computed for both conditions and segmented into conventional EEG bands for analysis ($\delta$, $\theta$, $\alpha$, $\beta$, $\gamma$). The plot reveals a significant decrease in spectral power in the $\delta$, $\theta$, $\alpha$ bands under the K+ conductance reduction condition, indicated by the triple asterisks (t-test, *p* > 0.001). Shaded areas represent the standard deviation across simulations and channels, while solid lines represent the average PSD. Right) These heatmaps illustrate the normalized differences in spectral power for the $\delta$, $\theta$, $\alpha$ and $\beta$ bands, demonstrating the robustness of the observed effect in the cortical-reduced Lausanne463 connectome. No significant difference was found in the gamma bands in any of the tested configurations. The x- and y-axes represent variations in the inhibitory and excitatory ($E_{l_i}$ and $E_{l_e}$) reversal potentials, respectively. The highlighted cells correspond to the parameters used in the figure on the left. The color gradient indicates the magnitude and direction of the normalized spectral power difference, with red signifying a positive difference and blue a negative difference (Student's t-tests, **: *p* < 0.001, **: *p* < 0.01, *: *p* < 0.05, ns: *p* > 0.05). This figure demonstrates that the decrease in the $\delta$, $\theta$, $\alpha$ channels is a robust finding that holds true across a range of reversal potential values for both excitatory and inhibitory neural populations using the cortical-reduced Lausanne463 connectome.</figcaption></figure>