Abstract
Psychoactive drugs can transiently perturb brain physiology while preserving brain structure. The role of physiological state in shaping neural function can therefore be investigated through neuroimaging of pharmacologically induced effects. Previously, using pharmacological neuroimaging, we found that neural and experiential effects of lysergic acid diethylamide (LSD) are attributable to agonism of the serotonin2A receptor (Preller et al., 2018). Here, we integrate brainwide transcriptomics with biophysically based circuit modeling to simulate acute neuromodulatory effects of LSD on human cortical largescale spatiotemporal dynamics. Our model captures the interareal topography of LSDinduced changes in cortical blood oxygen leveldependent (BOLD) functional connectivity. These findings suggest that serotonin2Amediated modulation of pyramidalneuronal gain is a circuit mechanism through which LSD alters cortical functional topography. Individualsubject model fitting captures patterns of individual neural differences in pharmacological response related to altered states of consciousness. This work establishes a framework for linking molecularlevel manipulations to systemslevel functional alterations, with implications for precision medicine.
Introduction
What are the respective roles of structure and physiology in shaping the spatiotemporal dynamics of networked neural systems? Areal differences in gene transcription, cellular architecture, and longrange connectivity patterns have been linked to areal differences in specialization of cortical function (Felleman and Van Essen, 1991; Hilgetag et al., 2016; Burt et al., 2018; Huntenburg et al., 2018; Demirtaş et al., 2019). Yet it remains unclear how embedded neurophysiological changes impact functional network organization and structurefunction relationships. Noninvasive neuroimaging of pharmacologically induced changes in brain function permits selective, transient, and in vivo perturbation of physiology, while preserving brain structure such as interareal axonal projections. Pharmacological neuroimaging therefore provides a tractable and valuable testbed for investigating the roles of neurophysiological properties in shaping patterns of largescale brain function.
The central role of physiological state in shaping functional brain dynamics and influencing human consciousness and behavior is supported by neuroimaging studies of the acute functional disturbances induced by psychopharmacological compounds. Serotonergic hallucinogens, in particular, produce rapid and profound alterations of consciousness that are linked to widespread, stereotyped patterns of functional network disruption (Tagliazucchi et al., 2014; Tagliazucchi et al., 2016; Lord et al., 2018; Muthukumaraswamy et al., 2013; Atasoy et al., 2017; Preller et al., 2018; Preller et al., 2020). Recently, our group has shown that lysergic acid diethylamide (LSD)induced disruptions of blood oxygen leveldependent (BOLD) functional connectivity (FC) and concomitant changes in consciousness are extinguished by pretreatment with ketanserin, a selective serotonin2A (5HT_{2A}) receptor antagonist (Preller et al., 2018). This result suggests that the 5HT_{2A} receptor plays a critical role in LSD’s mechanism of action (Vollenweider and Preller, 2020; Holze et al., 2021). Preller et al., 2018 analyzed regional changes in mean BOLD FC, called global brain connectivity (GBC), which is a graphtheoretic statistic that can be interpreted as a measure of functional integration (Cole et al., 2010; Anticevic et al., 2014). We found that the regional topography of LSDinduced changes in GBC is aligned with the topography of regional expression levels of HTR2A, the gene which encodes the 5HT_{2A} receptor (Preller et al., 2018). However, the circuit mechanism through which LSD alters cortical GBC topography remains unclear.
One approach to bridge this mechanistic gap is to develop biophysically based models of largescale brain activity that incorporate key features of neuronal and synaptic dynamics (Breakspear, 2017; Deco et al., 2011). These models, which are grounded in human neurobiology, can be first calibrated to healthystate data and then systematically perturbed through biophysically interpretable parameters. In doing so, selective manipulations at the synapse can be linked to their manifestations at empirically resolvable scales (Shine et al., 2021) – for instance, at the length scale and timescales probed by functional magnetic resonance imaging (fMRI). Moreover, recent advances in highthroughput transcriptomics have led to rich genomewide atlases of gene expression levels mapped across the human brain (Hawrylycz et al., 2012; Hawrylycz et al., 2015). Insofar as protein levels increase with the number of gene transcripts in a region (Liu et al., 2016), gene expression maps provide a principled way to infer the spatial distribution of proteins of interest, for example, drug targets. Gene expression maps therefore provide a means to simulate the regionally heterogeneous impacts of a drug on local circuit properties (Murray et al., 2018).
We extended this computational modeling approach to generate mechanistic insight into the topography of LSDinduced GBC alterations (Preller et al., 2018). GBC is pharmacologically elevated in sensory cortex, especially in visual cortex, and reduced in association cortex (Preller et al., 2018). Moreover, GBC changes correlate significantly with changes in subjects’ experience of consciousness, as determined by validated psychometric instruments (Studerus et al., 2010; Preller et al., 2018). Here, we investigated potential circuit mechanisms underlying these effects by extending a previously validated model of largescale neural dynamics (Deco et al., 2013; Deco et al., 2014; Demirtaş et al., 2019). We found that our model accurately captures the spatial topography of LSDinduced GBC changes. Our findings suggest that the distribution of 5HT_{2A} receptors is critical for generating the cortical topography of functional disruptions, and that neural gain is preferentially modulated on cortical pyramidal neurons, consistent with known neurobiology. Fitting the model at the individualsubject level, we found that the model captures patterns of individual differences in drug response that predict altered states of consciousness. Our work therefore advances pharmacological modeling from global effects at the group level to spatially heterogeneous effects in individuals. Broadly, this work establishes a flexible conceptual framework for linking mechanistic molecular hypotheses to human neuroimaging markers.
Results
To investigate the circuit mechanism through which LSD alters cortical GBC topography, we extended a previously validated biophysically based model of largescale neural circuit dynamics (Wong and Wang, 2006; Deco et al., 2013; Demirtaş et al., 2019; Figure 1). The model comprises many local microcircuits, each consisting of coupled excitatory and inhibitory neuronal populations. Population dynamics are driven by recurrent synaptic interactions that are governed by neurophysiologically interpretable parameters. Our model consists of 180 interconnected nodes, each of which represents one lefthemispheric cortical parcel in the Human Connectome Project’s (HCP’s) MultiModal Parcellation (MMP1.0) (Glasser et al., 2016). Note that we model only one cortical hemisphere, as the pharmacological fMRI data informing the model are approximately bilaterally symmetric (Preller et al., 2018), and the gene expression data are sampled from the left hemisphere but show bilateral correspondence (Hawrylycz et al., 2015). Nodes in the network interact through structured longrange excitatory projections. The relative strengths of these projections are informed by a diffusion MRIderived, groupaveraged (N = 339) structural connectivity (SC) matrix (Figure 1A; Van Essen et al., 2013; Demirtaş et al., 2019). The SC matrix is rowwise normalized such that total longrange inputs to each node are balanced and independent of parcel size. To evaluate model outputs against fMRIderived data, simulated synaptic activity in each node is transformed into an observable BOLD signal using the BalloonWindkessel model for the hemodynamic response (Friston et al., 2003; Stephan et al., 2007; Deco et al., 2014; Demirtaş et al., 2019; Figure 1B). We henceforth refer to this hemodynamically coupled neural circuit model as the ‘baseline’ or ‘unperturbed’ model.
To simulate the action of LSD in human cortex, we systematically perturb the baseline model. In earlier empirical analyses, we found that LSD’s influence on brain function and on behavior is primarily mediated through its agonist activity at the 5HT_{2A} receptor (Preller et al., 2018). Stimulation of 5HT receptors alters the response properties and membrane excitability of mammalian cortical neurons, such that firing rates are selectively enhanced for strongly depolarizing currents – that is, such that neural gain is enhanced (Araneda and Andrade, 1991; Andrade, 2011; Zhang and Arsenault, 2005). This 5HT receptormediated enhancement or modulation of gain is thought to be mediated specifically by the 5HT_{2A} receptor (Zhang and Arsenault, 2005). We therefore simulate the action of LSD by manipulating neural gain equations in the model (Figure 1C). Neural gain in the model is defined by a nonlinear expression of the form $r=f(aI)$ (Equation 2), which relates presynaptic current $I$ to postsynaptic firing rate $r$ by the nonlinear function $f$ with scalar gain parameter $a$. We mathematically express gain modulation as a fractional change in $a$, that is, $a\to (1+\delta ){a}_{0}$, where a_{0} is the unperturbed parameter value and δ determines the magnitude of modulation.
To account for the heterogeneous spatial distribution of 5HT_{2A} receptor density across cortical areas, gain modulation in each node is scaled in proportion to the local expression level of HTR2A – the gene which encodes the 5HT_{2A} receptor – using transcriptomic data from the Allen Human Brain Atlas (AHBA) (Figure 1—figure supplement 1; Hawrylycz et al., 2012; Burt et al., 2018). In cortex, 5HT_{2A} receptors are present on glutamatergic projection neurons and GABAergic interneurons (Mengod et al., 2010; Burnet et al., 1995; Santana et al., 2004). We introduce an additional degree of freedom into the model such that gain modulation can differentially impact the excitatory and inhibitory neuronal populations, thereby not making a priori assumptions as to the celltype specificity of 5HT_{2A}mediated gain modulation (Figure 1C). Specifically, the gain parameter $a$ of neuronal population $p$ in brain region $i$ is given by ${a}_{i}^{p}=(1+{h}_{i}{\delta}^{p}){a}_{0}^{p}$, where h_{i} is proportional to the expression level of gene HTR2A in the $i$ th region. For brevity, in what follows we drop the subscript $i$ and superscript $p$ notation. Note that the scalar gain parameter $a$ depends linearly on δ, and the gain function $r=f(aI)$ scales nonlinearly with $a$; thus, the gain functions scale nonlinearly with δ.
Quantifying functional organization with GBC maps
Functional organization of cortical dynamics, in the model and in the data, can be operationalized by constructing cortical GBC maps, for comparison to our prior empirical findings (Preller et al., 2018). The change in GBC induced by LSD provides a topographic map across cortical regions of the mean effect on FC, which we have previously compared to the topography of the HTR2A gene expression map (Preller et al., 2018). The first step in computing a GBC map is to construct an FC matrix, defined as the matrix of Pearson correlation coefficients between pairs of brain regions’ BOLD signal time traces. We apply the Fisher rtoZ transformation to offdiagonal elements of this matrix such that values are approximately normally distributed. Offdiagonal elements are then averaged along either rows or columns (as the FC matrix is symmetric), yielding a map of GBC which consists of a single scalar value for each region.
As reported in Preller et al., 2018, GBC maps were constructed from restingstate BOLD fMRI scans for 24 study participants. Subjects served as their own controls. GBC maps were constructed for each subject in both drug conditions (placebo and LSD). The groupaveraged contrast map of the change in GBC for LSD vs. placebo drug conditions – henceforth referred to as the ΔGBC map – specified the target model output (Figure 1D). In other words, we sought to recapitulate in our model the pattern of functional alterations indicated by changes in cortical GBC topography.
Empirically, the groupaveraged ΔGBC map reveals widespread LSDinduced hyperconnectivity of sensorysomatomotor regions, and hypoconnectivity of association regions (Preller et al., 2018). In the models, we generate a GBC map for both the unperturbed and neuromodulated models, corresponding to models for the placebo and LSD conditions, respectively (described in more detail below). From these maps, we compute a model ΔGBC map which is then quantitatively compared to the empirical ΔGBC map (Figure 1D).
Model captures the spatial topography of LSDinduced changes in cortical GBC
The operating point of the unperturbed model was first calibrated such that it best captured empirical FC in the placebo condition. This was achieved by performing a onedimensional grid search over the value of a global longrange coupling parameter, which uniformly scales the strengths of all longrange connections between nodes (Deco et al., 2014). We determined the value of this coupling parameter that yielded the greatest Spearman rank correlation between the offdiagonal elements of the simulated and empirical FC matrices. The coupling parameter was subsequently fixed at this value. The model FC matrix was then used to compute a baseline model GBC map – our computational analog of the empirical GBC map in the placebo condition.
Next we performed a twodimensional grid search over the two free model parameters which set the strength of gain modulation on excitatory and inhibitory neurons (${\delta}^{E}$ and ${\delta}^{I}$), such that a modulation of gain by 1% corresponds to $\delta =0.01$ (Figure 2A). For each combination of parameters (i.e., at each point on the parameter space grid), we updated the simulated gain functions and recomputed the system’s stable fixed point. This nowperturbed system is used to generate a new FC matrix, from which we again compute a model GBC map – this time, our computational analog of the empirical GBC map in the LSD condition. For each combination of free parameters, we then computed a model ΔGBC map by calculating the difference between the perturbed and unperturbed model GBC maps. To evaluate the quality of the model for each set of parameters, we evaluated the statistical similarity of each model ΔGBC map to the empirical ΔGBC map. We quantified similarity by computing the ‘loading’ of the model ΔGBC map onto the empirical ΔGBC map, which we define as the Cartesian dot product between the two ΔGBC maps (expressed as vectors), divided by the squared norm of the empirical ΔGBC map (the natural scale in the quantity). This metric was chosen specifically for its sensitivity to both the topography and magnitude of the model perturbation.
We found that the model ΔGBC map loaded most strongly onto the empirical ΔGBC map across a quasidegenerate range of parameter values (Figure 2A, offdiagonal yellow band). The location of this regime suggests that neural gain is preferentially modulated on excitatory pyramidal neurons, which nonetheless remain in a specific ratio with inhibitory interneurons. To investigate this further, we next characterized the relationship between changes in model excitatorytoinhibitory (E/I) balance and modelempirical loading. For each combination of gain modulation parameters (i.e., at each location on the twodimensional grid in parameter space), we computed the model’s E/I firing rate ratio, relative to its unperturbed value, where we define E/I ratio as the ratio of the mean excitatory firing rate (computed across nodes) to the mean inhibitory firing rate. We find that modelempirical loadings, when plotted as a function of fractional change in E/I ratio, collapse to form an approximately onedimensional curve which peaks at positive shifts in E/I ratio (Figure 2B). This onedimensional collapse suggests that gain modulationinduced changes in E/I balance are the primary mechanism through which GBC topography is altered under LSD. Moreover, the fact that the peak occurs for positive shifts in E/I balance suggests that LSD has a netexcitatory effect on cortical activity.
The model generates a strong functional perturbation that is aligned with the data, but does it also capture finegrained functional networkspecific and parcellevel effects? Empirically, LSD exhibits bidirectional effects across sensory and association brain networks, inducing hyperconnectivity of constituent brain regions in sensory and somatomotor networks, while inducing hypoconnectivity of regions in associative networks (Preller et al., 2018). To investigate network effects in the model, we set the gain modulation parameters to the values that yielded the greatest modelempirical loading (indicated by the black star in Figure 2A). We then characterized the distribution of model ΔGBC values across parcels within each functional network defined in the ColeAnticevic Brain Network Parcellation (CABNP) (Figure 2C; Ji et al., 2019). We found that empirical network specificity of GBC changes was recapitulated in the model (Figure 2D–E). This finding indicates that the strong modelempirical loading is not driven by select functional networks or by a small subset of all cortical parcels. In fact, at the level of parcels we find a remarkably strong spatial correspondence between ΔGBC topographies in the model and data (${r}_{s}=0.73$; p < 10^{−5}; Spearman rank correlation) (Figure 2F). Note that brain map values, here and in subsequent figures, are plotted on a flattened (i.e., unfolded) representation of the cortical surface, to better illustrate maps’ contiguous spatial topographies. This statistical association was significantly stronger than the association between the HTR2A gene expression map and the empirical ΔGBC map (p = 0.011; test for difference between dependent correlations) (Figure 2H). Heterogeneous physiological responses to LSD thereby play a critical role in shaping the cortical topography of LSDinduced GBC changes.
Global signal regression (GSR) remains a debated preprocessing step that is frequently included in fMRI analyses (Murphy and Fox, 2017). Preller et al., 2018 extensively characterized these data, with and without GSR applied. There, a statistically significant brainbehavior relationship between the neural and experiential effects of LSD was observed when – and only when – GSR was performed. Moreover, GSR was found to substantially alter the topography of the empirical ΔGBC map: when GSR was not used, the spatial correspondence of the ΔGBC map with the HTR2A expression map is dramatically attenuated (Preller et al., 2018). Here, we find that when GSR is not performed, maximal modelempirical loading and Spearman rank correlation between the model and empirical ΔGBC maps are both reduced significantly (Figure 2—figure supplement 1). This finding further supports the notion that GSR attenuates nonneural components or artifacts which otherwise obscure meaningful neuronal effects in pharmacological fMRI signals.
Receptor topography plays a critical role in shaping largescale functional effects
In addition to stimulating 5HT_{2A} receptors, LSD exerts predominately agonistic activity at (5HT)2C, 1A/B, 6, and 7 receptors, as well as dopamine D1 and D2 receptors (Nichols, 2004; Passie et al., 2008). Molecular signaling across these disparate receptor subtypes would impact neuronal physiology differently; however, these neuromodulatory systems are known to be important regulators of cortical gain (Thurley et al., 2008; Ferguson and Cardin, 2020). We therefore asked in our model whether modulation of gain via these other receptor subtypes also explains the spatial topography of cortical GBC changes. To test this, we repeated our grid search over model parameters, except that instead of modulating gain in proportion to HTR2A expression levels, we modulated gain in proportion to the expression level of, in turn: serotonergic genes HTR1A, HTR2C, and HTR7, and dopaminergic genes DRD1 and DRD2 (excluding 5HT_{1B} and 5HT_{6} encoding genes for lack of reliable transcriptomic data in the AHBA). The result of this process is illustrated in Figure 3A. Of the receptorencoding genes that we tested, using the HTR2A map in the model yields the maximal modelempirical loading and Spearman rank correlation between empirical and simulated ΔGBC maps (Figure 3B). These findings highlight the importance of the spatial distribution of drug targets in mediating their largescale functional effects.
We examined whether a strong modelempirical loading could be simply explained by general statistical properties of the HTR2A map, rather than its specific topographic pattern. To quantify the significance of the HTR2A map’s topography, we sought to establish statistical expectations for modelempirical loading under an appropriate null hypothesis. A widely adopted approach to quantifying significance for a complex, modelgenerated statistic (such as modelempirical loading) is to perform a nonparametric permutation test: by randomly permuting the brain map and regenerating a quantity of interest, one can construct a null distribution of expected outcomes due to random chance (Deco et al., 2018). When permutation tests of this type are applied to brain maps, the implicit null hypothesis is that any brain map with the same distribution of values but with a random topography is statistically likely to have produced a comparable or more extreme effect.
One important problem with this approach is that spatial autocorrelation – a characteristic statistical property of brain maps – is necessarily destroyed by random spatial permutations. Spatial autocorrelation is fundamentally important for two reasons: (i) it violates the assumption that samples are independent or exchangeable, an assumption which underlies many common statistical tests (including the permutation test); and (ii) this violation dramatically inflates pvalues in studies of brain maps (Burt et al., 2020). In light of this limitation, here we leverage a recent method to generate autocorrelationpreserving surrogate brain maps (Burt et al., 2020). By construction, these surrogate brain maps have randomized topographies but a fixed spatial autocorrelation (Figure 3C). Building a null distribution from these autocorrelationpreserving brain maps facilitates a test that controls for this important property.
To perform such a test, we repeatedly redefined gain functions in the model, each time by modulating gain in proportion to regional values of a random spatial autocorrelationpreserving surrogate brain map. Each surrogate map was constructed to have the same spatial autocorrelation structure as the HTR2A map. For each surrogate map, we regenerated a model ΔGBC map, from which we computed a modelempirical loading. The resulting null distribution of modelempirical loadings was then used to evaluate the significance of the HTR2A map’s specific topography (Figure 3D). We found that modulating gain in proportion to HTR2A expression levels produces a statistically significant loading (p = 0.008). The null distribution further reveals that none of the alternative gene expression maps reach significance of p < 0.05. These findings further implicate the 5HT_{2A} receptor in LSD’s mechanism of action, while simultaneously acting as a statistical control for the model complexity.
Model captures experientially relevant modes of neural variation
Our model captures cortical GBC alterations at the group level, but is it expressive enough to capture individual differences? To characterize individual variation – first, in the empirical data – we performed a principal components analysis (PCA) on subjects’ individual ΔGBC maps. We focused on the leading principal component (PC1) of ΔGBC variation (Figure 4A). By definition, PC1 corresponds to the spatial map that (linearly) captures maximal variance in ΔGBC across subjects. We found that PC1 correlates strongly with the groupaveraged ΔGBC map (${r}_{s}=0.72$; p < 10^{−5}; Spearman rank correlation). This indicates that individual differences are primarily driven by the strength with which subjects exhibit the groupaveraged pattern.
To explore individual variation in the model, for each subject, we repeated the twodimensional grid search over the two gainmodulating parameters. Modelempirical loading for each subject was computed using that subject’s ΔGBC map as the target model output. For each subject, we selected the combination of gain modulation parameters that maximized modelempirical loading. We then used those parameters to construct subjects’ model ΔGBC maps (Figure 4—figure supplement 1). We performed PCA on these subjectspecific maps to derive PC1 of ΔGBC variance in the model (Figure 4B). We found that, as in the empirical data, PC1 is topographically aligned with the grouplevel ΔGBC map (${r}_{s}=0.87$; p < 10^{−5}; Spearman rank correlation). Moreover, model PC1 and empirical PC1 exhibited similar spatial topographies (${r}_{s}=0.63$; p < 10^{−5}; Spearman rank correlation). These findings reveal a convergence between the dominant modes of neural variation in the data and in the model.
The model PC variance spectrum shows that variation across subjects’ simulated ΔGBC maps is almost entirely captured by the two leading model PCs (Figure 4C). That model variation lies within a twodimensional subspace is expected due to our fitting of two free parameters to ΔGBC. We also find topographic alignment between model PC2 and empirical PC2 (${r}_{s}=0.51$; p < 10^{−5}; Spearman rank correlation) (Figure 4D). However, the empirical PC variance spectrum reveals that the data exhibit approximately five significant modes of variation (Figure 4C). The data therefore exhibit three more modes than the model. A linear subspace analysis revealed that 67% of model variance fell within the subspace defined by the data’s leading five PCs (see Materials and methods). In other words, we found that the majority of model variation lies within a lowdimensional subspace defined by the principal modes of empirical variation.
With only two significant modes of variation, is the model expressive enough to capture individual differences in LSDinduced alterations of consciousness? As reported by Preller et al., 2018, changes in conscious experience were determined psychometrically using the fivedimensional altered states of consciousness (5DASC) questionnaire (Dittrich et al., 2010), which was administered to study participants 180 min after drug infusion. The short version of the 5DASC questionnaire quantifies changes in conscious experience along five distinct experiential dimensions: experiences of disembodiment, elemental imagery, changed meaning of percepts, blissful state, and spiritual experience – each of which is profoundly altered under LSD (Figure 4E). We performed a neurobehavioral linear regression analysis to find patterns of altered GBC that were predictive of experiential effects. Specifically, for each (parcel, experiential dimension) pair, we performed a linear regression across subjects, such that each sample in the regression between ΔGBC in parcel $i$ (the predictor variable) and change along experiential dimension $j$ (the target variable) corresponds to a (ΔGBC_{i}, ΔScore_{j}) pair for one subject.
For each experiential dimension, we performed one regression per parcel and aggregated the linear regression coefficients across all parcels to create a spatial brain map (Figure 4F). We henceforth refer to these maps as ‘experiential regression maps’. By construction, experiential regression maps reflect ΔGBC patterns that predict individuals’ change in experience. For instance, if a subject’s ΔGBC map strongly resembles the Meaning regression map (Figure 4F), then we expect that subject to have reported a relatively strong change in the meaning of percepts. How consistent are these experiential regression maps across experiential dimensions? We found evidence for a single dominant experiential regression map pattern (Figure 4G).
Returning to our original question – do modes of model variation have experiential relevance? – we compared PC1 maps (which characterize individual variation) to experiential regression maps (which are defined by their experiential relevance) (Figure 4G). We find notable correlations of model PC1 with the Meaning (${r}_{s}=0.38$; p < 10^{−5}; Spearman rank correlation) and Disembodiment (${r}_{s}=0.61$; p < 10^{−5}; Spearman rank correlation) experiential regression maps. Most strikingly, however, we find that the model PC1 map strongly resembles the dominant experiential regression map pattern (${r}_{s}=0.65$; p < 10^{−5}; Spearman rank correlation), which is defined as the first PC of the experiential regression maps. This indicates that the model captures patterns of functional individual variation that are linked to individual differences in pharmacologically induced changes in experience.
Discussion
In this study we integrated transcriptomic mapping into a biophysically based model of largescale cortical dynamics to investigate the circuit mechanisms underlying LSDinduced changes in cortical GBC topography. Recently, our group has shown that LSD induces experientially relevant changes in GBC through its agonist activity at the 5HT_{2A} receptor (Preller et al., 2018). Here, we found that the ΔGBC topography was captured in silico when LSD’s effect on molecular signaling was modeled as a modulation of neuronal gain, mediated by 5HT_{2A} receptor agonism, with preferential impact on pyramidal neurons. Moreover, our findings were specifically attributable to the topography of the HTR2A expression map – our transcriptomic proxy measure of the spatial distribution of 5HT_{2A} receptors. Finally, we showed that the model has enough expressivity to fit individual subjects and exhibits patterns of variation that predict druginduced alterations of consciousness. This work establishes a framework for linking molecularlevel manipulations to salient changes in brain function by integrating transcriptomics, biophysical modeling, and pharmacological neuroimaging.
Our findings complement a recent modeling study by Deco et al., 2018, which showed that a gainmodulatory mechanism captures changes in global spatiotemporal properties of grouplevel neural dynamics under LSD but did not examine cortical topographic effects as studied here. The authors of that study characterized changes in the distribution of pairwise correlations between slidingwindow dynamic FC matrices, or FC dynamics (FCD) (Hansen et al., 2015). FCD and other global statistical metrics have revealed interesting dynamical consequences of gain modulation in simulated neural systems (Shine et al., 2018; Deco et al., 2018; Li et al., 2019; Pfeffer et al., 2020). However, changes in global metrics may be driven by spatially nonspecific effects and by nonneuronal components or artifacts, particularly in imaging studies of serotonergic psychedelics (Lewis et al., 2017; Preller et al., 2018; Vollenweider and Preller, 2020). Our findings here demonstrate that gain modulation is a molecular mechanism that explains the specific spatial topography of LSD’s functionally disruptive effects. This key result illustrates how our proposed modeling framework could potentially be used in future studies to inform the development of therapeutics that precisely target pathological brain circuits while minimizing offtarget effects.
Whereas prior modeling work explored the effects of modulating neural gain uniformly across cell types (Deco et al., 2018), our study investigated the celltype specificity of LSD’s effects. Our finding that the best model fits are associated with elevation of E/I ratio and gain modulation that is biased toward excitatory populations suggests that LSD preferentially targets pyramidal neurons. This finding converges with multiple lines of experimental evidence in humans (Pazos et al., 1987; Hall et al., 2000), monkeys (Jakab and GoldmanRakic, 1998), and rats (Willins et al., 1997; Santana et al., 2004), which show that 5HT_{2A} receptors are preferentially found on apical dendrites of cortical pyramidal neurons. Moreover, 5HT_{2A} receptor agonists have been shown to induce substantial increases in cortical pyramidal neuron firing rates (MartínRuiz et al., 2001; Puig et al., 2003; LladóPelfort et al., 2018). Stateoftheart singlecell RNA sequencing data across multiple cortical areas in humans also reveals that HTR2A expression is significantly elevated in excitatory cell types, relative to inhibitory cell types (Figure 1—figure supplement 1).
Most prior models of psychedelics’ pharmacological effects in humans have been fitted to grouplevel statistical measures (Deco et al., 2018; Preller et al., 2019; Kringelbach et al., 2020; Herzog et al., 2020; Pfeffer et al., 2020). Here, we have shown that by fitting to individual subjects, the model can capture experientially relevant modes of neural variation. Our model fits the majority of subjects well (Figure 4—figure supplement 1). However, the number of subjects included in the empirical data collection supporting this study was not statistically well powered for individual differences predictive analyses, and several subjects were too idiosyncratic in their neural responses for us to accurately predict individuals’ experiential responses to LSD. Pharmacological modeling studies informed by other imaging modalities and powered by larger sample sizes are needed to fully characterize the predictive power and limitations of this approach.
LSD primarily mediates its effects on the human brain through 5HT_{2A} receptors (Preller et al., 2018; Holze et al., 2021), but LSD also acts as an agonist at other serotonergic and dopaminergic receptor sites (Nichols, 2004; Passie et al., 2008). In particular, LSD exhibits high affinity for 5HT_{1A} receptor subtypes, which function as inhibitory autoreceptors in the raphe nuclei and induce postsynaptic membrane hyperpolarization in limbic areas (Filip and Bader, 2009). Thus, 5HT_{2A} and 5HT_{1A} receptors mediate opposing influences on membrane excitability. Interestingly, this bidirectionality appears consistent with the ‘antipsychedelic’ properties of partial 5HT_{1A} receptor agonists in humans (Pokorny et al., 2016). In addition, 5HT_{2A} receptors are found in all cortical laminae (Jakab and GoldmanRakic, 1998; Burnet et al., 1995), often colocalized with 5HT_{1A} receptors on pyramidal neurons (Araneda and Andrade, 1991; AmargósBosch et al., 2004). 5HT_{2A} receptors are also found on GABAergic interneurons, though relatively infrequently (Willins et al., 1997; Santana et al., 2004).
Functional neuroimagingderived measures including GBC are sensitive to the global signal (GS), which represents shared variation across brain regions and may partially reflect nonneuronal physiological, movement and scannerrelated artifacts (Murphy et al., 2013; Murphy and Fox, 2017). In particular, GS artifacts exhibit dramatic differences in clinical populations and following pharmacological manipulations (Driesen et al., 2013; Yang et al., 2014; Power et al., 2017; Yang et al., 2017; Lewis et al., 2017). The data processing decision to include or forego GSR therefore has impact on findings in studies of psychedelics (Preller et al., 2018; Preller et al., 2020): it has been shown that results frequently do not replicate across samples when GSR is not performed (Tagliazucchi et al., 2016; Preller et al., 2018; Müller et al., 2017; Preller et al., 2020), while studies that use GSR have reported replicable findings (Preller et al., 2018; Preller et al., 2020). Here, we found that multiple measures of modelempirical similarity were significantly improved when GSR was performed (Figure 2—figure supplement 1), providing further support for its use. However, use of GSR remains debated (Fox et al., 2009; Saad et al., 2012), and in general there is no single ‘right’ way to process restingstate data (Murphy and Fox, 2017). In principle, it may be possible to circumvent some complications introduced by GSR by using a combination of spatial and temporal ICAbased denoising (Glasser et al., 2018), but this has not yet been tested in pharmacological fMRI studies.
In this study we leveraged recent advances in generative null modeling to establish statistical benchmarks while controlling for the confounding influence of spatial autocorrelation in analyses of brain maps (Burt et al., 2020). Spatial autocorrelation occurs ubiquitously in empirical brain maps and can dramatically inflate pvalues in analyses of both cortical and subcortical brain maps (Burt et al., 2020). Importantly, however, maps of brain features often exhibit the most marked differences across (as opposed to within) neuroanatomical structures. For instance, although gene expression profiles exhibit characteristic hierarchical patterns of intracortical variation (Burt et al., 2018), corticocortical gene expression variance is small relative to corticosubcortical variance (Hawrylycz et al., 2015). In addition to spatial autocorrelation within brain structures, these sharp distinctions between brain structures (such as cortex vs. thalamus) introduce additional bias into spatially naive permutation tests performed at the wholebrain level (e.g., as used by Deco et al., 2018).
We used gene expression maps, derived from the AHBA, as proxy correlates of receptor densities (Hawrylycz et al., 2012; Hawrylycz et al., 2015; Burt et al., 2018). Carlyle et al., 2017 found that variation in gene expression levels across human brain regions generally well captured regional variation in protein expression levels. Positron emission tomography (PET) is a neuroimaging modality which can measure regional variation in receptor availability in vivo. Beliveau et al., 2017 directly compared PET maps for serotonin receptors and found good alignment with gene expression levels from the AHBA. A key limitation of using PETderived maps for model simulations of pharmacology is that for many receptors of interest, PET mapping is not yet possible due to lack of developed radioligands (Pike, 2016). Our study demonstrates proof of concept that transcriptomic mapping can inform largescale models of pharmacological neuroimaging.
Our parsimonious computational model of LSD’s effects in cortex treats excitatory and inhibitory neurons as statistical ensembles, providing a description of neural dynamics at a level of resolution and complexity appropriate for proofofconcept comparisons with BOLD neuroimaging data. However, future work which expands, constrains, and tests biophysical models of pharmacology will be critical to address a number of open questions that remain. Neuromodulation and disease processes may predominately influence spatiotemporal properties within rather than across areas, suggesting the need for multiscale models that go beyond descriptions of brain areas in aggregate. Because our model is defined at the level of cortical parcels, it cannot speak to changes in connectivity that occur over smaller spatial scales, particularly among neurons within the parcels themselves. Our findings indicate that this coarse dynamical description is sufficient to capture regional GBC differences, but future work that goes beyond regional meanfield modeling may be needed to fully resolve the finegrained effects of pharmacology on within and betweenregion FC. Celltype specific pharmacological effects can also be further investigated by multicompartment models or models which include multiple interneuron subtypes. Celltype specific effects could be introduced into such models by integrating bulk transcriptomics with singlecell RNAsequencing to inform the joint distribution of pharmacological targets across distinct areas and cell types, respectively (Lake et al., 2016; Burt et al., 2018). Moreover, while 5HT_{2A} receptors are predominately found in cortex (Figure 1—figure supplement 1A–C), future model extensions with subcortical structures can incorporate the 5HT_{2A} receptorrich claustrum (Hall et al., 2000), which may play a role in mediating the effects of serotonergic psychedelic drugs (Barrett et al., 2020). Comparing pharmacological models to data from multiple imaging modalities may also generate new and complementary insights: for instance, the mechanisms underlying LSDinduced changes in broadband oscillatory power (Muthukumaraswamy et al., 2013; CarhartHarris et al., 2016) and increased neural signal diversity (Schartner et al., 2017) remain unclear.
Recent years have experienced a resurgence of clinical interest in the use of psychedelics as therapeutics for the treatment of mood disorders, alcohol and substance abuse, and endoflife distress in terminallyill patients (CarhartHarris and Goodwin, 2017; Nichols et al., 2017; Vollenweider and Preller, 2020). The common therapeutic mechanism of serotonergic psychedelic drugs remains unclear, although multiple lines of evidence now point to 5HT_{2A} receptormediated glutamate release (Vollenweider and Kometer, 2010; Mason et al., 2020; Vollenweider and Preller, 2020). Spatial gradients in drug targets, such as serotonergic receptor subtypes, can be approximated using the topography of gene transcripts – particularly in the case of receptors for which the field lacks suitable PET radioligands. Studies that leverage transcriptomic mapping (Grandjean et al., 2021) and regionally heterogeneous modeling (Demirtaş et al., 2019) will further elucidate the mechanisms underlying the complex functional signatures of pharmacologically induced neuromodulatory effects (Pfeffer et al., 2020; Alamia et al., 2020).
Materials and methods
Structural neuroimaging data
Request a detailed protocolGroupaveraged lefthemispheric SC matrices were constructed from diffusion MRI data using probabilistic tractography for 339 unrelated subjects from the HCP 900subject data release (Van Essen et al., 2013). SC matrices were parcellated into 180 areas using the HCP’s MMP1.0 (Glasser et al., 2016). Diagonal elements of the SC matrix were set identically to zero, as the dynamical model (described below) explicitly includes selfcoupling terms. Moreover, the SC matrix was rowwise normalized such that the total longrange inputs to each node were normalized. This normalization procedure instantiates the assumption that each local microcircuit receives a balance of local and longrange inputs. The CABNP was used to determine parcels’ functional network assignments (Ji et al., 2019). Pairwise interparcel geodesic distance matrices, which were used for constructing spatial autocorrelationpreserving surrogate maps (Burt et al., 2020), were computed by averaging over surfacebased distances between grayordinate vertices in each pair of parcels $i$ and $j$, where surfacebased distances were computed across the lefthemisphere midthickness surface in the HCP atlas.
Pharmacological fMRI data
Request a detailed protocolA complete description of data collection and preprocessing procedures was reported in Preller et al., 2018. Briefly, fMRI data derived from a doubleblind, placebocontrolled, and withinsubject study design, ensuring that each subject served as their own control. Twentyfive healthy human participants received either (i) placebo or (ii) LSD (100 µg po). Restingstate scans were collected 75 min following drug administration. fMRI data were obtained and processed following HCPcompliant acquisition standards (Glasser et al., 2013). Behavioral measures were derived from the short version of the 5DASC questionnaire (Dittrich et al., 2010), which includes 45 items that comprise five different scales: spiritual experience, blissful state, disembodiment, elementary imagery, and changed meaning of percepts. The 5DASC questionnaire was administered 180 min following drug infusion. One subject was excluded from analysis due to failure in registration caused by an improper head position. Thus, 24 study participants were analyzed in this study.
Transcriptomic data
Request a detailed protocolThe AHBA is a publicly available transcriptional atlas of microarray data containing samples from hundreds of neuroanatomical structures in six normal postmortem human brains. Microarray expression data and all accompanying metadata were downloaded from the AHBA (http://human.brainmap.org) (Hawrylycz et al., 2012; Hawrylycz et al., 2015). An exhaustive description of our preprocessing procedure was reported in Burt et al., 2018. Briefly, cortical microarray data in volumetric space were mapped onto subjects’ native twodimensional cortical surfaces by minimizing threedimensional Euclidean distance between microarray samples and grayordinate vertex coordinates. For parcels that were not directly sampled, we performed a surfacebased interpolation of sample expression values using a Voronoi diagram approach to create a dense gene expression map (at the level of grayordinates); this interpolated map was then parcellated by averaging across grayordinate vertices in each parcel. One representative microarray probe was selected for each unique gene, and gene expression profiles for selected gene probes were zscored. Finally, groupaveraged expression profiles were computed for genes whose expression profiles were reliable in at least four of the six subjects. These steps yielded groupaveraged gene expression values across 180 lefthemispheric cortical areas.
Singlenucleus transcriptomic data that quantify the expression of HTR2A across distinct cell types in human cortex was obtained from the Allen Brain Institute’s Human Multiple Cortical Areas SMARTseq database (https://portal.brainmap.org/atlasesanddata/rnaseq/humanmultiplecorticalareassmartseq). We downloaded the CSV file containing ‘Gene expression by cluster, trimmed means’ and aggregated columns that included either ‘Exc’ or ‘Inh’.
Modeling largescale neural dynamics
Request a detailed protocolWe adapted the biophysically based largescale computational model described in Deco et al., 2014. This model reduces the complexity and the number of local microcircuit parameters in a spiking neural network model using a dynamical meanfield approach (Wong and Wang, 2006). By leveraging a statistical description of ensemble neural activity and exploiting the long time constants of NMDA receptors, the model reduces a highdimensional spiking neural network to a computationally tractable twodimensional dynamical system. Each node in the model comprises recurrently coupled excitatory (E) and inhibitory (I) populations, the dynamics of which are described via coupled nonlinear stochastic differential equations. Timevarying activity of the excitatory and inhibitory synaptic currents in brain region $i$, denoted respectively by ${I}_{i}^{E}$ and ${I}_{i}^{I}$, is defined as
where ${I}_{b}$ is constant background input current; $S$ is the synaptic gating variable, that is, the fraction of open channels; $G$ is a global conductance parameter which scales the strengths of longrange connections; $J$ is the effective NMDA conductance; $C$ is the SC matrix; and remaining parameters are constants. Synaptic current ${I}^{p}$ for population $p\in \{E,I\}$ is transformed into firing rate ${r}^{p}$ via the transfer function
where $a$, $b$, and $d$ are celltype specific parameters governing the form of the inputoutput relation. In particular, note that parameter $a$, which regulates the slope of the transfer function, corresponds to the neural gain parameter. Note that while this function is unbounded and therefore does not saturate, we confirmed that the nodeaveraged firing rates in the gainmodulated model (which do not exceed ∼15 Hz throughout the parameter sweep in Figure 2A) remain in a neurobiologically plausible firing rate regime where this approximation of the FI curve does not break down.
Finally, the time evolution for synaptic gating variable $S$ is given by
where τ is a synaptic time constant, γ is a kinetic parameter, and σ is the standard deviation of the stochastic Gaussian input noise ν (Wong and Wang, 2006; Deco et al., 2014). Following Deco et al., 2014, simulated synaptic covariance matrices were approximated analytically by linearizing Equations 1 to 3 around the dynamical system’s stable fixed point.
Neural gain modulation was implemented in the model by introducing a parametric rescaling of gain parameter $a$ for each cell type, in proportion to regional gene expression levels of HTR2A:
where h_{i} denotes the expression level of HTR2A in brain region $i$, ${\delta}^{E}$ and ${\delta}^{I}$ regulate the strength of the perturbation applied to excitatory and inhibitory cell populations, respectively. Thus, twodimensional parameter sweeps reported in this study, for instance in Figure 2A, correspond to sweeps over free parameters ${\delta}^{E}$ and ${\delta}^{I}$. Where brain maps (gene expression maps as well as surrogate maps) were used to modulate gain, they were first linearized using the error function (following Demirtaş et al., 2019).
Following Deco et al., 2014, we implemented feedback inhibition control (FIC) in the model such that the firing rates ${r}_{i}^{E}$ were constrained to be 3 Hz prior to neural gain modulation. This is accomplished by tuning the weight ${w}^{EI}$, which regulates the strength of local inhibitorytoexcitatory feedback, within each simulated brain region. Crucially, we note that FIC was calculated prior to but not following gain modulation, under the assumption that FIC represents neurobiological selfregulatory processes which balance neural activity in the brain on timescales much longer than the characteristic timescales of transient pharmacological manipulations.
Modeling the hemodynamic response
Request a detailed protocolTo facilitate model comparisons with empirical fMRI data, excitatory synaptic activity within each brain region (given by ${S}^{E}$ in Equation 3) was transformed to a BOLD signal using the BalloonWindkessel model (Friston et al., 2003). In this model, the hemodynamic response obeys the following coupled system of equations:
where $x$ denotes the vasodilatory signal, $f$ the blood inflow, $v$ the blood volume, and $q$ the deoxyhemoglobin content; and parameters ρ, τ, κ, γ, and α are the resting oxygen extraction fraction, hemodynamic transit time, rate of signal decay, rate of flowdependent elimination, and the Grubb’s exponent, respectively (Friston et al., 2003). In turn, the corresponding BOLD signal $y$ is calculated as
where V_{0} is the resting blood volume fraction. The three dimensionless magnetic field strengthdependent parameter values k_{1}, k_{2}, and k_{3} were derived for a magnetic field strength of 3 T using Appendix A of Heinzle et al., 2016; all other hemodynamic parameter values were taken from Obata et al., 2004. Simulated BOLD covariance matrices were derived by linearizing these equations and then algebraically transforming the linearized covariance matrix, using a procedure which we previously reported in Demirtaş et al., 2019. Specifically, we semianalytically derive the covariance matrix $P$ of the dynamical system by numerically solving the Lyapunov equation:
where $A$ is the Jacobian matrix and ${Q}_{N}$ is the noise covariance matrix. Note that this expression is solved for the full sixdimensional dynamical system which includes two synaptic variables, ${S}^{E}$ and ${S}^{I}$, as well as four hemodynamic state variables $x$, $f$, $v$, and $q$. The linearized BOLD covariance matrix is then given by
where $K$ is a matrix of partial derivatives of the BOLD signal with respect to the six dynamical variables (Demirtaş et al., 2019).
Values for all static parameters in the model are provided in Table 1.
Modeling fitting
Request a detailed protocolTuning the model parameters consisted of two key sequential steps. First, the dynamical operating point of the model was calibrated such that (prior to gain modulation) the model optimally captured empirical FC derived from the placebocondition fMRI scan. To this end, we performed a onedimensional grid search over the global coupling parameter $G$, which linearly scales all longrange interactions between nodes (Equation 1). The grid search was over the range [0.01, 0.85] with a step size of 0.01. For each value, we solved for the stable fixed point of the system and then computed the linearized BOLD covariance matrix around that point. From the BOLD covariance matrix, we computed the FC matrix using the transformation from covariance to Pearson correlation coefficient. We computed the Spearman rank correlation between the $N=180(1801)/2$ uppertriangular elements of the model FC matrix and the placebocondition empirical FC matrix. At the group level (i.e., using the groupaveraged empirical FC), this step yielded an optimal global coupling value of $G=0.85$, which resulted in a Spearman rank correlation of ${r}_{s}=0.45$. This value of $G$ was used for all groupaveraged analyses in this study. We note that if the upper end of the range was extended beyond $G=0.85$, the maximal fit to the grouplevel data occurred at $G=0.89$. However, at this value the fit increases only marginally ($\delta {r}_{s}<{10}^{3}$), while the system is then positioned so close to a critical bifurcation that even modest modulations of gain will dynamically destabilize the system. This class of computational models is known to produce optimal fits near the edge of stability (see, e.g., Deco et al., 2013).
After calibrating the system, next we performed a twodimensional grid search over model parameters $\delta E$ and $\delta I$, which determine the strength of excitatory and inhibitory gain modulation, respectively. For each combination of parameters, we updated the neural gain equations (Equation 5) in the model. We solved for the new fixed point and computed the linearized BOLD covariance matrix around that point. FC and GBC were then computed as described above. The difference between model GBC maps (expressed as Fisher Zvalues) – perturbed minus baseline – yielded a model ΔGBC map. We then computed the modelempirical loading between the model ΔGBC map and the empirical ΔGBC map (LSD minus placebo). As illustrated in the heatmaps, the grid search for the gain modulation parameters was over the range [0, 0.03] with a step size of 0.003. Gain parameters that maximized modelempirical loading were then used for all analyses involving a model ΔGBC map, either at the group or subject level.
Global brain connectivity
Request a detailed protocolEmpirical GBC was computed using inhouse Matlab tools for all grayordinates in the brain using an HCPharmonized version of the FreeSurfer software (Fischl et al., 2002; Cole et al., 2011; Anticevic et al., 2013; Anticevic et al., 2014), as reported in Preller et al., 2018. Briefly, for each subject, we computed the fMRI timeseries correlation between pairs of grayordinates, thereby constructing an FC matrix. Offdiagonal FC matrix elements were Fisher Ztransformed and then averaged across rows to compute GBC values. Maps were then parcellated using the HCP’s MMP1.0 (Glasser et al., 2016). This yielded a parcellated GBC map for each subject where each parcel’s value represents the mean FC of that region to all other regions. Parcellated GBC maps were averaged across subjects to construct the groupaveraged GBC map. Similarly, to compute model GBC maps, we used the model to generate a parcellated BOLD FC matrix. Offdiagonal FC matrix elements were Fisher Ztransformed and then averaged across rows to compute GBC values. All illustrated GBC map values (empirical and simulated) are Fisher Zvalues.
Global signal regression
Request a detailed protocolGSR was performed on empirical fMRIderived timeseries data using widely adopted approaches (Anticevic et al., 2013; Cole et al., 2011), as described in Preller et al., 2018. GSR was also performed in the model, as GSR not only removes artifactual signal components but neuronal signal as well. In the model, however, we performed GSR directly on the analytically approximated BOLD covariance matrices using the following approach. Let y_{i} be the BOLD signal in brain region $i$, and let $x$ be the GS – that is, the mean graymatter BOLD signal across regions. Note that $x$ and $y$ are both implicit functions of time. We construct a linear regression model of the form
where the residual signal following GSR is given by
What we seek is an expression for the covariance between GSregressed BOLD signal residuals in regions $i$ and $j$, that is,
where $\overline{\u03f5}$ denotes the temporal average of $\u03f5$. Substituting Equations 12 and 13, this expression can be rewritten as
By manipulating Equation 12, β can be expressed in terms of covariances:
Substituting this expression into Equation 15 and simplifying, we obtain
From the definition of GS $x={\sum}_{i}{y}_{i}$, Equation 17 can be rewritten as
Collecting the summations and noting that
we obtain
Equation 20 analytically relates the covariance matrix of the raw BOLD signal (righthand side) to the covariance matrix of GSregressed residuals (lefthand side). We used this expression to perform GSR on model BOLD covariance matrices. FC matrices (i.e., Pearson correlation matrices) were then computed from GSregressed BOLD covariance matrices.
Generative modeling of surrogate brain maps
Request a detailed protocolSpatial autocorrelationpreserving surrogate maps were constructed following Burt et al., 2020 using the Pythonbased BrainSMASH toolbox (https://brainsmash.readthedocs.io/). This approach operationalizes spatial autocorrelation in a brain map by computing a variogram, which quantifies variance in the brain map as a function of pairwise distance between areas. To construct each surrogate map, the empirical brain map is iteratively shuffled (to randomize topography), smoothed (to reintroduce spatial autocorrelation), and rescaled (to recover empirical spatial autocorrelation). After performing this procedure, we resampled our surrogate map values from the empirical brain map such that the distribution of values does not change: each surrogate map can then be conceptualized as one random realization of an (approximately) spatial autocorrelationpreserving permutation of the original brain map.
For each surrogate map $\overrightarrow{h}$, the surrogate map values $\{{h}_{i}\}$ were substituted into Equation 5 such that the modulation of gain in each region was scaled in proportion to the values of the surrogate map, rather than the empirical HTR2A map. This ‘surrogate model’ was then used to recompute the statistic under consideration. Repeating for and aggregating across all surrogate maps yields a null distribution for the expected value of the statistic under the null hypothesis that the outcome of the statistical test is an artifact of spatial autocorrelation, and therefore not unique to a specific spatial topography. Specifically, where we report pvalues that derive from a spatial autocorrelationpreserving surrogate map test, we are reporting the proportion of samples in the null distribution which were more extreme than the test statistic.
Experiential regression maps
Request a detailed protocolTo identify characteristic patterns of changes in GBC that predict subjectlevel changes in 5DASC scores, we constructed spatial brain maps of linear regression coefficients (i.e., beta coefficient maps or simply ‘beta maps’). To calculate the beta map value in the $i$ th cortical parcel for experiential dimension $d$, we solved for ${\beta}_{i}^{d}$ in the linear regression defined by
where $\delta {s}^{d}$ is the change in score $\delta s$ along experiential dimension $d$; α and β are constant and linear regression coefficients, respectively; and $\mathrm{\Delta}{\mathrm{GBC}}_{i}$ is the change in GBC for parcel $i$; and vectors include 24 elements (one per subject).
Linear decomposition
Request a detailed protocolTo define lowdimensional linear subspaces within the 180dimensional neural state space, we used PCA to determine the dominant spatial modes of variation across different collections of brain maps. We let $\mathbf{\mathbf{X}}$ be the $M\times N$ matrix comprised of $M$ brain maps, each represented as an $N$dimensional vector, where columns of $\mathbf{\mathbf{X}}$ have been meansubtracted. The spatial covariance matrix is then given by the $N\times N$ symmetric matrix $\mathbf{\mathbf{C}}$, where
The symmetric matrix $\mathbf{\mathbf{C}}$ can in general be decomposed according to
where $\mathbf{\mathbf{P}}$ is an $N\times N$ orthogonal matrix whose $i$ th column $\overrightarrow{{p}_{i}}$ is an eigenvector of $\mathbf{\mathbf{C}}$, and $\mathbf{\mathbf{\Lambda}}$ is an $N\times N$ diagonal matrix whose $i$ th diagonal element ${\lambda}_{i}$ is the eigenvalue associated with eigenvector p_{i}. Eigenvector $\overrightarrow{{p}_{i}}$ is a PC of $\mathbf{\mathbf{X}}$, and the normalized eigenvalue ${\lambda}_{i}/\mathrm{Tr}(\mathbf{\mathbf{\Lambda}})$ is the fraction of total variance in $\mathbf{\mathbf{X}}$ that occurs along $\overrightarrow{{p}_{i}}$. We will assume without loss of generality that eigenvectors are ordered such that ${\lambda}_{i+1}\le {\lambda}_{i}$. In other words, we assume that eigenvectors are arranged such that eigenvalues are placed in descending order. Finally, note that if $M<N$, then $NM+1$ eigenvalues are identically zero (i.e., $\mathbf{\mathbf{C}}$ will not be full rank).
To determine which PCs were statistically significant, we performed simple permutation testing. Each row in $\mathbf{\mathbf{X}}$ (i.e., each brain map) was randomly permuted to obtain a new matrix ${\mathbf{\mathbf{X}}}^{\prime}$. The procedure described above was then performed on ${\mathbf{\mathbf{X}}}^{\prime}$ and the fraction of variance captured per PC (i.e., $\overrightarrow{{\lambda}^{\prime}}/\mathrm{Tr}({\mathbf{\mathbf{\Lambda}}}^{\prime})$) was recorded. To construct the gray dashed lines in Figure 4, we performed this procedure $N=1000$ times and then computed the fifth percentile of the resulting null distribution (for each individual PC).
To investigate shared dimensions of variability in the model and in the data, we first performed PCA on both the model and empirical subjectspecific ΔGBC maps. Permutation testing revealed that the leading two model PCs and the leading five empirical PCs were statistically significant. We then derived the covariance matrix $\mathbf{\mathbf{C}}$ for the $24\times 180$ data matrix $\mathbf{\mathbf{X}}$, comprised of subjects’ model ΔGBC maps, per Equation 22. Total model variance $V$ was determined by computing the trace of $\mathbf{\mathbf{C}}$. Let $\mathbf{\mathbf{P}}$ denote the $180\times 5$ matrix whose columns correspond to the leading five empirical PCs; then we performed a linear transformation of $\mathbf{\mathbf{C}}$ into the fivedimensional subspace defined by the columns of $\mathbf{\mathbf{P}}$ according to
Finally, we determined the total variance ${V}^{\prime}$ contained in the subspace by computing the trace of ${\mathbf{\mathbf{C}}}^{\prime}$. Then the fraction of model variance that falls within the data’s fivedimensional subspace (67%) is given by ${V}^{\prime}/V$.
Data and code availability
Request a detailed protocolA Pythonbased implementation of the model is provided as an openaccess software package, BRAINTRIPS (BRainwide Activity Induced by Neuromodulation via TRanscriptomicsInformed Pharmacological Simulation): https://github.com/murraylab/braintrips (Burt, 2021; copy archived at swh:1:rev:2f4c9bb63ec01c40e4374e277d03c7ff28e92713). All processed neuroimaging data needed to reproduce main findings in this study are made publicly available on BALSA (https://balsa.wustl.edu/study/show/5XMxv).
Data availability
A Pythonbased implementation of the model is provided as an openaccess software package, BRAINTRIPS (BRainwide Activity Induced by Neuromodulation via TRanscriptomicsInformed Pharmacological Simulation): https://github.com/murraylab/braintrips (copy archived at https://archive.softwareheritage.org/swh:1:rev:2f4c9bb63ec01c40e4374e277d03c7ff28e92713). All processed neuroimaging data needed to reproduce main findings in this study are made publicly available on BALSA: https://balsa.wustl.edu/study/show/5XMxv.

BALSAID 5XMxv. Transcriptomicsinformed modeling of pharmacological neuroimaging effects of LSD.
References

A HighResolution in vivo atlas of the human brain's Serotonin SystemThe Journal of Neuroscience 37:120–128.https://doi.org/10.1523/JNEUROSCI.283016.2016

Dynamic models of largescale brain activityNature Neuroscience 20:340–352.https://doi.org/10.1038/nn.4497

The therapeutic potential of psychedelic drugs: past, present, and futureNeuropsychopharmacology 42:2105–2113.https://doi.org/10.1038/npp.2017.84

A multiregional proteomic survey of the postnatal human brainNature Neuroscience 20:1787–1795.https://doi.org/10.1038/s4159301700112

Emerging concepts for the dynamical organization of restingstate activity in the brainNature Reviews Neuroscience 12:43–56.https://doi.org/10.1038/nrn2961

How local excitationinhibition ratio impacts the whole brain dynamicsJournal of Neuroscience 34:7886–7898.https://doi.org/10.1523/JNEUROSCI.506813.2014

Book5dasc: Questionnaire for the assessment of altered states of consciousness. A short introductionZurich, Switzerland: PSIN PLUS.

The impact of NMDA receptor blockade on human working memoryrelated prefrontal function and connectivityNeuropsychopharmacology 38:2613–2622.https://doi.org/10.1038/npp.2013.170

Mechanisms underlying gain modulation in the cortexNature Reviews Neuroscience 21:80–92.https://doi.org/10.1038/s415830190253y

The global signal and observed anticorrelated resting state brain networksJournal of Neurophysiology 101:3270–3283.https://doi.org/10.1152/jn.90777.2008

Canonical genetic signatures of the adult human brainNature Neuroscience 18:1832–1844.https://doi.org/10.1038/nn.4171

LargeScale gradients in human cortical organizationTrends in Cognitive Sciences 22:21–31.https://doi.org/10.1016/j.tics.2017.11.002

Effects of hallucinogens on neuronal activityCurrent Topics in Behavioral Neurosciences 36:75–105.https://doi.org/10.1007/7854_2017_473

Handbook of Behavioral Neuroscience123–138, Distribution of 5ht receptors in the central nervous system, Handbook of Behavioral Neuroscience, Elsevier, 10.1016/S15697339(10)700746.

Increased thalamic restingstate connectivity as a core driver of LSDinduced hallucinationsActa Psychiatrica Scandinavica 136:648–657.https://doi.org/10.1111/acps.12818

Biophysical modeling of LargeScale brain dynamics and applications for computational PsychiatryBiological Psychiatry: Cognitive Neuroscience and Neuroimaging 3:777–787.https://doi.org/10.1016/j.bpsc.2018.07.004

Broadband cortical desynchronization underlies the human psychedelic stateJournal of Neuroscience 33:15171–15183.https://doi.org/10.1523/JNEUROSCI.206313.2013

HallucinogensPharmacology & Therapeutics 101:131–181.https://doi.org/10.1016/j.pharmthera.2003.11.002

Psychedelics as medicines: an emerging new paradigmClinical Pharmacology & Therapeutics 101:209–219.https://doi.org/10.1002/cpt.557

The pharmacology of lysergic acid diethylamide: a reviewCNS Neuroscience & Therapeutics 14:295–314.https://doi.org/10.1111/j.17555949.2008.00059.x

Considerations in the development of reversibly binding PET radioligands for brain imagingCurrent Medicinal Chemistry 23:1818–1869.https://doi.org/10.2174/0929867323666160418114826

Psilocybin induces TimeDependent changes in global functional connectivityBiological Psychiatry 88:197–207.https://doi.org/10.1016/j.biopsych.2019.12.027

Enhanced repertoire of brain dynamical states during the psychedelic experienceHuman Brain Mapping 35:5442–5456.https://doi.org/10.1002/hbm.22562

Dopamine increases the gain of the inputoutput response of rat prefrontal pyramidal neuronsJournal of Neurophysiology 99:2985–2997.https://doi.org/10.1152/jn.01098.2007

The neurobiology of psychedelic drugs: implications for the treatment of mood disordersNature Reviews Neuroscience 11:642–651.https://doi.org/10.1038/nrn2884

Psychedelic drugs: neurobiology and potential for treatment of psychiatric disordersNature Reviews Neuroscience 21:611–624.https://doi.org/10.1038/s4158302003672

A recurrent network mechanism of time integration in perceptual decisionsJournal of Neuroscience 26:1314–1328.https://doi.org/10.1523/JNEUROSCI.373305.2006

Altered global signal topography in schizophreniaCerebral Cortex 27:5156–5169.https://doi.org/10.1093/cercor/bhw297

Gain modulation by serotonin in pyramidal neurones of the rat prefrontal cortexThe Journal of Physiology 566:379–394.https://doi.org/10.1113/jphysiol.2005.086066
Decision letter

Timothy E BehrensSenior and Reviewing Editor; University of Oxford, United Kingdom

James M ShineReviewer; University of Sydney, Australia

John GriffithsReviewer
Our editorial process produces two outputs: i) public reviews designed to be posted alongside the preprint for the benefit of readers; ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.
Acceptance summary:
This paper will be of interest to scientists working on computational modelling of neuroimaging data, and on the neural effects of psychedelic drugs and other pharmacological interventions. The study is wellmotivated. The statistical and data analytic methodologies are rigorous and advanced. The with conclusions are wellsupported by the presented data. The modelling methodology includes technical innovations that are potentially of broad utility and importance.
Decision letter after peer review:
Thank you for submitting your article "Transcriptomicsinformed largescale cortical model captures topography of pharmacological neuroimaging effects of LSD" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Timothy Behrens as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: James M Shine (Reviewer #1); John Griffiths (Reviewer #2).
The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.
Essential revisions:
None. We like the manuscript and would be happy to publish it tomorrow if you choose. However, we have made some suggestions below that we think would improve the manuscript, which you can implement at your discretion.
Reviewer #1 (Recommendations for the authors):
P4 – perhaps better to state that GBC 'can be interpreted' as a measure of functional integration.
P4 – I realise it is passe to ask for citations of one's own work, however the authors may wish to cite a recent review on the use of neural mass models linking structure and function (Shine et al., 2021 in Nature Neuroscience), as the topic of the review is stronglyaligned with the authors approach.
P6 – is it realistic to model gain as a continually increasing function? There is a natural ceiling to how high the firing rate of a neuron can be, which suggests that a sigmoid transfer function might be a more sensible function. The authors could mitigate this concern by confirming that the firing rate of their neural populations is bounded within reasonable limits by other features of their model (e.g., EI balance).
P7 – are the offdiagonal elements of the FC matrix normalized separately for 5HT2A vs. placebo conditions? If so, did the authors first check to determine whether there were systematic differences between 5HT2A and placebo that may have been diminished through normalisation?
P8 – the authors conclude that the model fits reflect the fact that "neural gain is preferentially modulated on excitatory pyramidal neurons", however I wonder whether a more parsimonious description would be that "neural gain is preferentially modulated on excitatory pyramidal neurons, which nonetheless remain in a specific ratio with inhibitory interneurons". Note that this result is further expanded in Figure 2B, but I worry that the interim conclusion may mislead from the final result.
P8 – I really liked the utilisation of other 5HT and DA receptor maps and permutation testing.
Reviewer #2 (Recommendations for the authors):
It would be informative to include some discussion of the connectivity normalization choice. The connectivity matrix diagonals are set to zero and each row is rescaled to unity, which is equivalent to using the Laplacian. The result of this is that every brain region receives an identical total level of input from other regions in the network. This is particularly interesting given the principal metric of interest is (changes in) global brain connectivity (row/column averages of the FC matrix). The equivalent maps for the anatomical connectivity will be uniform, for the reasons detailed above. Could the authors please discuss: what, if any, is the neurobiological, and/or mathematical rationale for this normalization choice, and what are their thoughts on the above considerations.
Reviewer #3 (Recommendations for the authors):
I really enjoyed reading this manuscript, it is very wellwritten and easy to follow. The study appears to be methodologically sound, however, given my lack of direct expertise in the modeling of dynamical systems, I am will refrain from giving more specific comments and/or suggestions for improving the scientific quality of these analyses.
https://doi.org/10.7554/eLife.69320.sa1Author response
Reviewer #1 (Recommendations for the authors):
P4 – perhaps better to state that GBC 'can be interpreted' as a measure of functional integration.
We now state: “…global brain connectivity (GBC), which is a graphtheoretic statistic that can be interpreted as a measure of functional integration.”
P4 – I realise it is passe to ask for citations of one's own work, however the authors may wish to cite a recent review on the use of neural mass models linking structure and function (Shine et al., 2021 in Nature Neuroscience), as the topic of the review is stronglyaligned with the authors approach.
This new review is now cited.
P6 – is it realistic to model gain as a continually increasing function? There is a natural ceiling to how high the firing rate of a neuron can be, which suggests that a sigmoid transfer function might be a more sensible function. The authors could mitigate this concern by confirming that the firing rate of their neural populations is bounded within reasonable limits by other features of their model (e.g., EI balance).
We now state in the Methods: “Note that while this function is unbounded and therefore does not saturate, we confirmed that the nodeaveraged firing rates in the gainmodulated model (which do not exceed ∼15Hz throughout the parameter sweep in Figure 2A) remain in a neurobiologically plausible firingrate regime where this approximation of the FI curve does not break down.”
P7 – are the offdiagonal elements of the FC matrix normalized separately for 5HT2A vs. placebo conditions? If so, did the authors first check to determine whether there were systematic differences between 5HT2A and placebo that may have been diminished through normalisation?
We now include this clarification: “The location of this regime suggests that neural gain is preferentially modulated on excitatory pyramidal neurons, which nonetheless remain in a specific ratio with inhibitory interneurons.”
P8 – the authors conclude that the model fits reflect the fact that "neural gain is preferentially modulated on excitatory pyramidal neurons", however I wonder whether a more parsimonious description would be that "neural gain is preferentially modulated on excitatory pyramidal neurons, which nonetheless remain in a specific ratio with inhibitory interneurons". Note that this result is further expanded in Figure 2B, but I worry that the interim conclusion may mislead from the final result.
P8 – I really liked the utilisation of other 5HT and DA receptor maps and permutation testing.
Reviewer #2 (Recommendations for the authors):
It would be informative to include some discussion of the connectivity normalization choice. The connectivity matrix diagonals are set to zero and each row is rescaled to unity, which is equivalent to using the Laplacian. The result of this is that every brain region receives an identical total level of input from other regions in the network. This is particularly interesting given the principal metric of interest is (changes in) global brain connectivity (row/column averages of the FC matrix). The equivalent maps for the anatomical connectivity will be uniform, for the reasons detailed above. Could the authors please discuss: what, if any, is the neurobiological, and/or mathematical rationale for this normalization choice, and what are their thoughts on the above considerations.
We now state in the Methods: “Diagonal elements of the SC matrix were set identically to zero, as the dynamical model (described below) explicitly includes selfcoupling terms. Moreover, the SC matrix was rowwise normalized such that the total longrange inputs to each node were normalized. This normalization procedure instantiates the assumption that each local microcircuit receives a balance of local and longrange inputs.”
https://doi.org/10.7554/eLife.69320.sa2Article and author information
Author details
Funding
National Institute of Mental Health (R01MH112746)
 John D Murray
National Institute of Mental Health (DP5OD012109)
 Alan Anticevic
National Institute of Mental Health (R01MH108590)
 Alan Anticevic
Swiss National Science Foundation (P2ZHP1\161626)
 Katrin H Preller
Swiss Neuromatrix Foundation (20150103)
 Franz X Vollenweider
Usona Institute (20152056)
 Franz X Vollenweider
National Institute on Alcohol Abuse and Alcoholism (P50AA01287016)
 John H Krystal
 Alan Anticevic
Brain and Behavior Research Foundation (Independent Investigator Grant)
 Alan Anticevic
Simons Foundation (Pilot Award)
 Alan Anticevic
 John D Murray
National Center for Advancing Translational Sciences (UL1TR000142)
 Alan Anticevic
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
This research was supported by grants from the NIH (R01MH112746, JDM; DP5OD012109, AA; R01MH108590, AA); the Swiss National Science Foundation (P2ZHP1_161626, KHP); the Swiss Neuromatrix Foundation (2015–0103, FXV); the Usona Institute (2015–2056, FXV); the NIAAA (P50AA01287016, AA and JHK); the NARSAD Independent Investigator Grant (AA), a SFARI Pilot Award (JDM, AA); and the Yale CTSA grant (UL1TR000142 Pilot Award, AA). The funding sources had no involvement in the study design; nor in the collection, analysis, and interpretation of data; nor in the writing of the manuscript; nor in the decision to submit the manuscript for publication.
Ethics
Clinical trial registration ClinicalTrials.gov: NCT02451072.
Human subjects: This study did not involve any new data collection. Secondary data analysis was performed on the data set described in our prior publication (Preller et al., 2018, eLife). The Swiss Federal Office of Public Health, Bern, Switzerland, authorized the use of LSD in humans, and the study was approved by the Cantonal Ethics Committee of Zurich (KEKZH_No: 2014_0496). The study was registered at ClinicalTrials.gov (NCT02451072).
Senior and Reviewing Editor
 Timothy E Behrens, University of Oxford, United Kingdom
Reviewers
 James M Shine, University of Sydney, Australia
 John Griffiths
Publication history
 Received: April 13, 2021
 Accepted: June 8, 2021
 Version of Record published: July 27, 2021 (version 1)
Copyright
© 2021, Burt et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 792
 Page views

 92
 Downloads

 2
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.