Transcriptomics-informed large-scale cortical model captures topography of pharmacological neuroimaging effects of LSD

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 serotonin-2A receptor (Preller et al., 2018). Here, we integrate brain-wide transcriptomics with biophysically based circuit modeling to simulate acute neuromodulatory effects of LSD on human cortical large-scale spatiotemporal dynamics. Our model captures the inter-areal topography of LSD-induced changes in cortical blood oxygen level-dependent (BOLD) functional connectivity. These findings suggest that serotonin-2A-mediated modulation of pyramidal-neuronal gain is a circuit mechanism through which LSD alters cortical functional topography. Individual-subject model fitting captures patterns of individual neural differences in pharmacological response related to altered states of consciousness. This work establishes a framework for linking molecular-level manipulations to systems-level functional alterations, with implications for precision medicine.

Trial registration: ClinicalTrials.gov NCT02451072.

Keywords: computational model; functional connectivity; gene expression; human; neuroscience; pharmacological neuroimaging.

Figure 1.. Schematic overview of the biophysical modeling framework.

(A) Each node in the large-scale model represents a cortical microcircuit comprised of recurrently coupled excitatory (E) and inhibitory (I) neuronal populations. The model includes one node for each of the 180 left-hemispheric parcels in the Human Connectome Project’s Multi-Modal Parcellation (MMP1.0). Nodes interact through structured long-range excitatory projections, the strengths of which are constrained by a diffusion magnetic resonance imaging (MRI)-derived structural connectivity (SC) matrix. (B) Simulated synaptic activity in each node is transformed to a simulated blood oxygen level-dependent (BOLD) signal using the Balloon-Windkessel model of the hemodynamic response. (C) Lysergic acid diethylamide (LSD)’s effect on cortical microcircuitry is modeled as a modulation of neural gain due to serotonin-2A (5-HT_2A) receptor agonism by the LSD molecule. The degree to which neural gain is modulated within an area is scaled in proportion to the regional expression level of HTR2A, the gene which encodes the 5-HT_2A receptor protein. Gain curves of the excitatory and inhibitory neuronal populations are modulated independently, permitting cell-type specific effects. (D) Global brain connectivity (GBC), a graph-theoretic functional measure, is dramatically altered following LSD administration. The functional MRI (fMRI)-derived map of the change in GBC (ΔGBC) under LSD, relative to placebo, specifies the target model output. To simulate brain function in the LSD and placebo drug conditions, we simulate GBC maps with and without gain modulation, respectively. We compute the difference between the model GBC maps to construct a simulated ΔGBC map. Quantitative comparisons between empirical and model ΔGBC maps determine how well the model captures the topography of LSD-induced functional disruptions.

Figure 1—figure supplement 1.. In the human brain, HTR2A is predominately expressed in cortical pyramidal neurons.

In the human brain, HTR2A is predominately expressed in cortical pyramidal neurons. (A) Cortical topography of the HTR2A expression map. (B) Whole-brain topography of the HTR2A expression map. For subcortex, we use the 358 subcortical parcels in the Cole-Anticevic Brain Network Parcellation (CAB-NP) (Ji et al., 2019). Expression levels are linearly rescaled such that the minimum value is zero, and the cortical parcel-wise average is one. The large difference between expression levels in cortex and subcortex is much greater than the variance across parcels within the cortex. Note that the Allen Human Brain Atlas has unilateral sampling of gene expression in the left hemisphere (Hawrylycz et al., 2015), and therefore the map is made bilaterally symmetric at the parcel level for cortex and coordinate level for subcortex. Gene expression mapping follows the method of Burt et al., 2018. (C) HTR2A expression levels grouped by gross anatomical structure. Box plots mark the median and inner quartile ranges for expression levels across parcels within each anatomical structure, and whiskers indicate the 95% confidence interval. ‘Subcortex (aggregated)’ comprises parcel expression levels for all subcortical structures (i.e., all 358 subcortical parcels). Expression of HTR2A is significantly higher in cortex than in subcortex (${\displaystyle W=57}$; p < machine precision; Wilcoxon signed-rank test). (D) The distribution of HTR2A expression levels across excitatory (red) and inhibitory (blue) human cortical cell types. HTR2A is significantly more expressed in excitatory neurons than in inhibitory neurons (W = 141,943; p < machine precision; Wilcoxon signed-rank test).

Figure 2.. HTR2A-mediated excitatory gain modulation captures effects of lysergic acid diethylamide (LSD) on human cortical global brain connectivity (GBC) topography.

(A) Two-dimensional grid search over the two free model parameters. Parameters govern the gain modulation of inhibitory and excitatory neuronal populations in the model. Model-empirical loading – the quantity shown in the heatmap – is defined as the dot product between the empirical change in GBC (ΔGBC) map and a model ΔGBC map, normalized by the squared norm of the empirical ΔGBC map. Loading is maximized for the combination of parameters indicated by the black star. (B) At each point on the grid (i.e., for each combination of gain-modulatory parameters), we computed the excitatory-to-inhibitory (E/I) firing rate ratio, expressed in terms of its unperturbed value. This amounts to computing $r^{\prime }/r$, where $r^{\prime }$ is the E/I ratio in the model with gain modulation, and $r$ denotes the E/I ratio in the model without gain modulation. E/I ratio is defined in the model as the ratio of the mean excitatory firing rate (computed across nodes) to the mean inhibitory firing rate. (C) Functional network assignments for each cortical parcel are determined by the Cole-Anticevic Brain Network Parcellation (CAB-NP): ventral multi-modal (VMM), language (LAN), dorsal attention (DAN), posterior multi-modal (PMM), primary visual (VIS), secondary visual (VIS2), frontoparietal (FPN), cingulo-opercular (CON), default mode (DMN), orbito-affective (OAN), auditory (AUD), and somatomotor (SOM) networks. Network colors mirror (Ji et al., 2019). (D) Functional network-level comparisons between simulated (solid) and empirical (striped) z-scored ΔGBC map values. Box plots mark the median and inner quartile ranges for parcels in each network, and whiskers indicate the 95% confidence interval. (E) Distributions of z-scored ΔGBC map values across cortical parcels in primary sensory networks (AUD, VIS, VIS2, SOM; black line with no fill) and association networks (gray fill with no line). Endpoints of the horizontal black lines (top) correspond to the distributions’ means. Distributions significantly differ empirically (p < 10⁻⁴; spatial autocorrelation-preserving surrogate map test) and in the model (p = 0.02; spatial autocorrelation-preserving surrogate map test). (F) Spatial topographies of the dense and parcellated empirical ΔGBC maps; the strongest-loading model ΔGBC map; and the HTR2A gene expression map. Maps are portrayed on flattened (unfolded) representations of the cortical surface. (G) Scatter plot illustrating the parcel-wise relationship between the strongest-loading model ΔGBC map and the empirical ΔGBC map. (H) Spearman rank correlations between the empirical ΔGBC map and: (i) the strongest-loading model ΔGBC map; and (ii) the HTR2A expression map. Empirical ΔGBC topography is better explained by the dynamical model than by the HTR2A map (p < 0.05; test for dependent correlations).

Figure 2—figure supplement 1.. Effect of global signal regression (GSR) on model performance.

Effect of GSR on model performance. (A) The HTR2A gene expression map is significantly correlated with the empirical GS-regressed change in global brain connectivity (ΔGBC) map (${\displaystyle r_s=0.40}$; p < 10⁻⁵; Spearman rank correlation). In contrast, the empirical ΔGBC map constructed without GSR correlates only very weakly, and negatively, with the HTR2A map (${\displaystyle r_s=-0.15}$; p = 0.05). (B) Model-empirical loadings for the empirical ΔGBC map with no GSR (Figure 2A). Black star indicates model parameters which yielded maximal loading. (C) Model-empirical loading as a function of fractional change in model excitatory-to-inhibitory (E/I) ratio following gain modulation (related to Figure 2B). (D) The strongest loading model ΔGBC maps, with and without GSR. (E) Performing GSR significantly improves: the Spearman rank correlation between the HTR2A expression map and the empirical ΔGBC map (left; p < 10⁻⁵; test for difference between dependent correlations); maximal model-empirical loading (center; p = 0.03; spatial autocorrelation-preserving surrogate map test); and the Spearman rank correlation between the model and empirical ΔGBC maps (right; p < 10⁻⁵; test for difference between dependent correlations). Purple, without GSR; cyan, with GSR.

Figure 3.. The topography of lysergic acid diethylamide (LSD)-induced cortical global brain connectivity (GBC) changes is specifically attributable to the spatial distribution of 5-HT_2A receptors, as indexed by the HTR2A map.

(A) Two-dimensional grid searches over free model parameters. For each heatmap, gain modulation is scaled in proportion to regional expression levels of different serotonergic (HTR) and dopaminergic (DRD) receptor-encoding genes, each of which is agonized by LSD. Black stars indicate maximal model-empirical loadings for each heatmap. (B) Model-empirical loading (left axis) and Spearman rank correlation (right axis) are greatest when gain is modulated by regional expression levels of HTR2A. Model change in GBC (ΔGBC) maps used in this analysis were generated using the gain-modulatory parameters that maximized model-empirical loadings, as indicated on each heatmap. (C,D) Following Burt et al., 2020, we generate surrogate brain maps with randomized spatial topographies that, by construction, exhibit spatial autocorrelation that has been matched to that of the HTR2A map. These spatial autocorrelation-preserving surrogate brain maps are used to construct a null distribution of the expected magnitude of model-empirical loading under random chance. Each sample in the null distribution (gray; N = 1000) was constructed by modulating gain in proportion to the values in a random spatial autocorrelation-preserving surrogate brain map, then re-computing model-empirical loading. Colored lines correspond to the different receptor-encoding genes (as reported in panel B). ** p < 0.01.

Figure 4.. Model fits to individual subjects capture experientially relevant modes of neural variation.

(A–D) Characteristic patterns of inter-individual variability in cortical change in global brain connectivity (ΔGBC) topography. (A) The empirical group-averaged ΔGBC map (top), and the leading principal component (PC1) computed across 24 subjects’ empirical ΔGBC maps (bottom). (B) The model ΔGBC map fit to the group-averaged data (top), and the PC1 computed across 24 subjects’ model ΔGBC maps (bottom). (C) Empirical (green) and model (purple) PC variance spectra. The first five empirical PCs and the first two model PCs survived permutation testing (p < 0.05, 1000 permutations). (D) Empirical and model PC2 maps. (E–H) Linking individual differences in ΔGBC to individual differences in lysergic acid diethylamide (LSD)-induced alterations of consciousness. (E) Changes in subjects’ conscious experience under LSD relative to placebo as determined by the five-dimensional (5D) altered states of consciousness (5D-ASC) questionnaire: disembodiment, elementary imagery, changed meaning of percepts, blissful state, and spiritual experience. Box plots mark the median and inner quartile ranges for each scale, and whiskers indicate the 95% confidence interval. Positive values indicate an increase under lysergic acid diethylamide (LSD). (F) Experiential regression maps are constructed by performing linear regressions, across subjects, between changes in an experiential score (the target variable) and changes in GBC within a single parcel (the predictor variable). For each 5D-ASC scale, we performed 180 regressions – one per parcel – across 24 subjects. Brain maps illustrate the first-order regression coefficients. Experiential regression maps therefore illustrate patterns of GBC variation that predict experiential variation. (G) We performed principal components analysis (PCA) on the five experiential regression maps to derive experiential regression map principal components (PCs). Top: variance spectrum for the experiential regression map PCs. PC1 survives permutation testing (p < 0.05, 1000 permutations). Bottom: experiential regression map PC1, which captures 60% of variance across experiential regression maps. (H) Spearman rank correlations between experiential regression maps and neural PC1 maps.

Figure 4—figure supplement 1.. Empirical change in global brain connectivity (ΔGBC) maps for each subject.

Empirical ΔGBC maps for each subject. For each subject, we repeated the model-fitting procedure that was used for the group-level analysis: the global coupling parameter was first set to the value that maximized the Spearman rank correlation between off-diagonal elements of the model functional connectivity (FC) matrix and empirical subject-specific placebo FC matrix. Gain-modulatory parameters were then set to the values that maximized model-empirical loading (computed using the subject-specific empirical map).