Abstract
Dynamic reconfigurations of the functional connectome across different connectivity states are highly heritable, predictive of cognitive abilities, and linked to mental health. Despite their established heritability, the specific polymorphisms that shape connectome dynamics are largely unknown. Given the widespread regulatory impact of modulatory neurotransmitters on functional connectivity, we comprehensively investigated a large set of single nucleotide polymorphisms (SNPs) of their receptors, metabolic enzymes, and transporters in 674 healthy adult subjects (347 females) from the Human Connectome Project. Preregistered modulatory neurotransmitter SNPs and dynamic connectome features entered a Stability Selection procedure with resampling. We found that specific subsets of these SNPs explain individual differences in temporal phenotypes of fMRI-derived connectome dynamics for which we previously established heritability. Specifically, noradrenergic polymorphisms explained Fractional Occupancy, i.e., the proportion of time spent in each connectome state, and cholinergic polymorphisms explained Transition Probability, i.e., the probability to transition between state pairs, respectively. This work identifies specific genetic effects on connectome dynamics via the regulatory impact of modulatory neurotransmitter systems. Our observations highlight the potential of dynamic connectome features as endophenotypes for neurotransmitter-focused precision psychiatry.
Significance Statement
Understanding how genetic variations affect brain activity and connectivity can unlock new insights into cognitive abilities and mental health. This study reveals that specific genetic variations influence how long the brain stays in different connectivity states and how it transitions between these states. These genetic variations were found in two modulatory neurotransmitter systems: acetylcholine and noradrenaline. These findings suggest that brain connectivity patterns influenced by genetics could serve as markers for personalized psychiatric treatment, pushing the boundaries of precision psychiatry.
Introduction
The time-varying characteristics of the brain's functional connectivity architecture (the connectome) are specific to the individual, highly heritable, and predictive of individual cognitive abilities (Vidaurre et al., 2017; Barber et al., 2021; Jun et al., 2022). Such dynamic characteristics consist of temporal features which depict the connectome's sequence of state (spatial connectivity patterns), i.e., the trajectory of the connectome through the space of possible states. These dynamic connectome features are associated with cognitive functioning and implicated in several psychiatric and neurological disorders (Wang et al., 2020; Yang et al., 2021; Chen et al., 2023). Importantly, particular temporal features of connectome dynamics, such as the proportion of the total recording time spent in each connectome state (Fractional Occupancy) and the probability to transition between specific pairs of connectome states (Transition Probability), have been linked to behavioral performance (Vidaurre et al., 2017; Eichenbaum et al., 2020; Jun et al., 2022) and found to be heritable (Vidaurre et al., 2017; Jun et al., 2022). Specifically, our previous work has established substantial genetic effects (h2 ∼40%) and strong behavioral relevance of Fractional Occupancy and Transition Probability (Jun et al., 2022). This work further highlighted that such heritability is specific to the temporal features and does not encompass the spatial connectivity pattern and topology (specifically modularity) of the connectome states. Despite such strong support for genetic effects on the temporal characteristics of functional connectome dynamics, it is unknown which specific genes carry this effect. The goal of the current study is to identify genetic variants associated with these dynamic connectome phenotypes for which heritability has been established.
Recent theoretical advances and empirical evidence suggest that ascending neuromodulatory inputs play a pivotal role in driving such dynamic fluctuations in connectome states and cognitive functioning. This role may be due to the widespread volume transmission and long-lasting regulatory impact of modulatory neurotransmitters (Pfeffer et al., 2018; Shine, 2019; Shine et al., 2019; Totah et al., 2019). Unlike wiring transmissions, which involves rapid point-to-point signaling, neuromodulators such as acetylcholine and monoamines are commonly released nonsynaptically and instead diffusing through the extracellular space. Such diffusion allows a large number of cells to detect these neuromodulators through extrasynaptic receptors. Indeed, pharmacological studies have shown that the release of all core modulatory neurotransmitters—i.e., acetylcholine (ACh; Klaassens et al., 2017), noradrenaline (NAd; Shine et al., 2018a), dopamine (DA; Shafiei et al., 2019), and serotonin (5-HT; Klaassens et al., 2017; Deco et al., 2018; Lord et al., 2019; Kringelbach et al., 2020)—can drastically alter large-scale connectome reconfigurations across various behavioral states. Further, spatial distributions of neurotransmitter receptors are aligned with the structural and intrinsic connectivity networks (Picard et al., 2013; Kringelbach et al., 2020; Hansen et al., 2022), facilitating the flow of information and supporting cognitive processes (Sadaghiani et al., 2017).
Single nucleotide polymorphisms (SNPs), the most common form of genetic variation (Shastry, 2009), have been extensively studied for their role in shaping individual differences in brain structure and function (Elliott et al., 2018). As specific and measurable genetic variants, SNPs enable fine-grained analyses of genetic contributions to brain function. SNPs can modify a gene's function by either altering its expression or by modifying the composition of its product, i.e., the protein, through changes of amino acids or splice variants, to name the most common means. Therefore, we propose that SNPs in the genes, encoding receptors, metabolic enzymes, and transporters of the modulatory neurotransmitter systems, are likely contributors to the heritability of dynamic connectome state transitions and associated cognitive functions.
If empirically confirmed, such genetic impact of neurotransmitter systems on connectome dynamics and cognition would have implications for understanding mental health. Connectome dynamics are implicated in numerous psychiatric and neurological conditions (Preti et al., 2017) and explain individual differences in certain cognitive abilities across different diagnostic categories (Zhu et al., 2021). Regarding Fractional Occupancy and Transition Probability in particular, associations have been observed with major depressive disorder (Wang et al., 2020), schizophrenia (Yang et al., 2021), and subjective cognitive decline (Chen et al., 2023). Moreover, the patho-etiology of most psychiatric disorders involves abnormalities of neurotransmitter-related proteins, such as receptors dysfunction and neurotransmitter imbalance (Mandal et al., 2022). Consequently, most pharmacological interventions for psychiatric disorders target modulatory neurotransmitter systems (Li et al., 2016, 2022). Not surprisingly, specific SNPs in these systems may contribute to individual differences in symptoms and treatment response (Kautzky et al., 2015; Castro-Martínez et al., 2018). Therefore, it is likely that modulatory neurotransmitter systems impact mental health and cognitive functioning through dynamic reconfigurations of the functional connectome and that polymorphisms in the former shape the latter.
Despite such strong translational implications for putative genetic effects on functional connectome dynamics, previous genetics neuroimaging studies have not investigated the dynamic connectome. These studies have primarily focused on the association between genetic polymorphisms and structural connectivity (Elliott et al., 2018; Nemmi et al., 2018), static (time-averaged) functional connectivity (Markett et al., 2016; Elliott et al., 2018; Simola et al., 2022; Wang et al., 2022), and more broadly cortical volume (Nemmi et al., 2018) and network activations (Whelan et al., 2012; Sadaghiani et al., 2017; Macare et al., 2018; Baker et al., 2019). Further, most of the previous work has separately examined the impact of each individual genetic polymorphism of interest on brain imaging phenotypes. However, in a complex biological system, phenotypic expression typically results from the joint contribution of multiple genetic factors. Therefore, our study aimed to fill a critical knowledge gap by demonstrating how genetic polymorphisms of the modulatory neurotransmitter systems jointly shape the heritable phenotypes of the functional connectome dynamics, i.e., Fractional Occupancy and Transition Probability (Jun et al., 2022).
To this end, we examined a comprehensive list of SNPs from genes encoding ionotropic and metabotropic receptors, metabolic enzymes, and transporters for each modulatory neurotransmitter system (ACh, NAd, DA, 5-HT). The quality control and inclusion/exclusion criteria of the SNPs as well as definition of connectome and cognitive phenotypes followed our preregistered project (https://doi.org/10.17605/OSF.IO/VF2ZW). For each neurotransmitter system, we employed a fivefold cross-validation approach with a stable feature selection algorithm (Stability Selection; Meinshausen and Bühlmann, 2010; Shah and Samworth, 2013; Hofner et al., 2015) to explain each connectome phenotype. This approach aimed to identify the joint contribution of all SNPs in each neurotransmitter system to the connectome dynamics phenotypes. Further, in light of the established relationship between the above-described connectome dynamics phenotypes and cognitive abilities (Vidaurre et al., 2017; Eichenbaum et al., 2020; Jun et al., 2022), we examined whether the same set of SNPs that explained connectome dynamics could concomitantly explain variability in cognitive abilities.
Materials and Methods
Figure 1 is a schematic representation of the overall approach and analysis subsections. SNP selection and phenotype definition followed our preregistered project. Note that statistical analyses deviated from the preregistered project because we identified a testing approach more pertinent to the goal of the present study (Ayers and Cordell, 2010), which we were not aware of at the time of preregistration (compare its successful application in comparable genotype–phenotype association studies in Valls et al., 2019; Zou et al., 2020).
An overview of the analysis pipeline. Panels A and B follow the steps outlined in our preregistered project to determine genotypes and phenotypes entering the study. Specifically, the top section of Panel A describes the process of selecting the single nucleotide polymorphisms (SNPs), while the bottom section details the imputation and quality control steps carried out on the Human Connectome Project (HCP) genotype data. The left side of Panel B depicts the hidden Markov model (HMM) estimates of temporal phenotypes (i.e., Fractional Occupancy and Transition Probability) of functional connectome dynamics (Jun et al., 2022). The right side of Panel B illustrates the factorization of the various cognitive measures obtained from the HCP dataset. Note that as a foundation for the current study, we have previously established the heritability of the connectome and cognitive phenotypes depicted in Panel B using a twin study approach in the same data (Jun et al., 2022). In Panel C, we divided the data into five subsets for each connectome phenotype and neurotransmitter system. Within the training set (comprising 4/5th of the full data), we employed Stability Selection approach, which utilizes elastic net regularization and resampling techniques, to select a subset of SNPs and covariates (age, sex, head motion, and 10 principal components representing the European-ancestry population structure) that best explain each temporal phenotype. Subsequently, we used this subset in the “full model” ([including the selected SNPs, covariates, and random genetic effects (RGX)] and the “base model” (including the same subset of covariates and RGX, but not SNPs) to be fitted to the test set. This process was repeated five times. For each connectome phenotype and neurotransmitter system, we combined the residuals from the full and base models across the five test sets and used a linear mixed-effects model to assess model performance (full vs base). In Panel D, we refitted the full models that outperformed their respective base models in Panel C to be fitted to each of the Cognitive Factors in the same test sets. Subsequently, we combined the residuals of the “cognitive full” and “cognitive base” models across the five test sets and used a linear mixed-effects model to evaluate model performance. Refer to Extended Data Figure 1-1 for a visualization of the population stratification of all HCP subjects with genotype data, Extended Data Figure 1-2 for the list of 60 selected genes from the Neurotransmitter Study of the Allen Institute for Brain Science, Extended Data Figure 1-3 for the 155 SNPs included in the analysis after quality control, Extended Data Figure 1-4 for a detailed visualization of the four HMM-derived states, and Extended Data Figure 1-5 for information on the 14 cognitive measures used in the analysis.
Figure 1-1
3D plot of top 3 principal components representing the population stratification of all 1,141 subjects. Genetic ancestry of subjects was inferred using SNPweights software package (see Materials and Methods 2.1). Download Figure 1-1, TIF file.
Figure 1-2
The list of select 60 genes from the Neurotransmitter Study of the Allen Institute for Brain Science. Complete details of the genes characterized by in situ hybridization (ISH) for the Neurotransmitter Study, conducted by the Allen Institute for Brain Science, are accessible via the following link (https://help.brain-map.org/download/attachments/2818165/HBA_ISH_GeneList.pdf?version=2&modificationDate=1614977648535&api = v2; see Table 3 in the link). From the comprehensive gene list provided in the above link, our study specifically included genes associated with all core neuromodulatory systems, specifically the acetylcholine, dopamine, serotonin, and epinephrine/norepinephrine neurotransmitter systems. ACh; acetylcholinergic, NAd; noradrenergic, DA; dopaminergic, 5HT; serotonergic, GPCR; G-protein coupled receptor. Download Figure 1-2, DOCX file.
Figure 1-3
List of 155 SNPs included in the analysis. Note that the variant type information of each variant (or SNP) has been acquired from (https://www.ensembl.org; accessed on 24 May 2023). CHR; chromosome number, ALT; alternate allele, REF; reference allele, ACh; acetylcholinergic, NAd; noradrenergic, DA; dopaminergic, 5HT; serotonergic. Download Figure 1-3, DOCX file.
Figure 1-4
HMM-derived connectome states. [A] HMM was used to estimate the connectome states that are common to all subjects investigated in our prior study (Jun et al., 2022). The z-scored functional connectivity (FC) matrix of each state is visualized. The rows and columns represent 139 regions organized according to their membership to canonical intrinsic connectivity networks (ICNs) (Yeo et al., 2011) including default mode network (DMN), frontoparietal network (FPN), dorsal attention network (DAN), salience/cingulo-opercular network (CON), sensory-motor network (SMN), visual network, and limbic network (Limbic; subcortical limbic regions omitted for visual succinctness). [B] A binary matrix showing a data-driven set of clusters of connections whose FC differed significantly across the four states. F values from a connection-wise ANCOVA of FC strength were threshold at 5% connection density followed by Network-Based Statistics (NBS) to control for multiple comparisons (Zalesky et al., 2010). Significant FC changes across the connectome states prominently include connections within the visual network, between DMN and numerous other ICNs, and between DAN and the visual network. Download Figure 1-4, TIF file.
Figure 1-5
List of 14 cognitive measures from HCP and their basic descriptions. Download Figure 1-5, DOCX file.
Subjects
We used genetics, resting-state fMRI and behavioral data from the Washington University-University of Minnesota (WU-Minn) consortium's Human Connectome Project (HCP) S1200 release (Van Essen et al., 2013). Participants were recruited, and informed consent was acquired by the WU-Minn HCP consortium according to procedures approved by the Washington University IRB (Glasser et al., 2013).
The genotyping data for 1,141 subjects was made available through the dbGAP repository. Briefly, the DNA samples were collected from either whole blood or saliva sample and genotyped using the Illumina Multi-Ethnic Global Array (MEGA) SNP array. This array included chip-specific content from PsychChip and ImmunoChip and provided comprehensive coverage of European, East Asian, and South Asian populations. To maintain ethnic homogeneity, only subjects with European genetic ancestry were included in the study using the genetic ancestry predicted with the SNPweights software package (Chen et al., 2013; Altmann and Mourao-Miranda, 2019; Extended Data Fig. 1-1). Consequently, 674 healthy adult subjects (aged 22–36 years, 347 females) with four complete resting-state fMRI scans (4,800 total timepoints) were included. Further details on the HCP data collection protocol (Van Essen et al., 2013; Smith et al., 2013a), cognitive measures (Barch et al., 2013), and inclusion and exclusion criteria (Van Essen et al., 2013; Smith et al., 2013a) can be found elsewhere.
Genotype imputation
Genotype imputation for the SNPs of interest (see below) was performed using minimac4 method and the HRC reference panel (version r1.1, which consists predominantly of European ancestry) through the Michigan Imputation Server (accessed on 29 Jan. 2022; Das et al., 2016). Prior to the imputation process, we conducted quality control of the HCP data using the toolbox provided by Will Rayner (http://www.well.ox.ac.uk/∼wrayner/tools/; accessed on Jan. 29, 2022). This step was taken to ensure the data was properly organized to be imputed with the HRC reference panel. The input data was prepared following the guidelines provided by the Michigan Imputation Server. The genotype imputation was done for each chromosome after several quality control and phasing steps by the server.
SNP selection and quality control
We targeted a specific set of SNPs in the genes that encode receptors, metabolic enzymes, ion channels, and transporters of modulatory neurotransmitters (i.e., ACh, NAd, DA, and 5-HT) in the human brain following our preregistered selection procedures. Specifically, we selected 60 genes related to these neurotransmitter systems from the Neurotransmitter Study (Extended Data Fig. 1-2) conducted by the Allen Institute for Brain Science. These genes were characterized by in situ hybridization (ISH) in the adult human brain. The SNPs associated with these genes were identified using the Variant Effect Predictor (VEP; https://useast.ensembl.org/info/docs/tools/vep/) tool. In particular, we identified SNPs with the upstream and downstream distance of 5,000 bp from the gene as a buffer (which is a default setting of VEP and a common choice in the field) meeting specific criteria, including having a Minor Allele Frequency (MAF) greater than 0.05, as well as being of either intron, synonymous, or missense variant types. These criteria were chosen to ensure that the selected SNPs represent common genetic variations and are pertinent to gene expression regulation and protein function. Further, we used GTEx Portal (https://gtexportal.org/; accessed on 02 Nov. 2021) to quantify each SNP's expression quantitative trait loci (eQTL) in cortical tissue, then selected SNPs with eQTL p value <0.05. This additional step allowed us to hone in on SNPs with substantial influence on gene expression in the cortex.
After identifying the SNPs to be included in the study, quality control measures were taken on the imputed HCP data using PLINK v2.0 (Purcell et al., 2007). SNPs with a MAF < 5%, a Hardy–Weinberg equilibrium (p ˂ 1 × 10−6), and missing call rates (>10%) across individuals were removed. Additionally, individuals with missing call rates (>20%) across the genotypes were excluded. Linkage disequilibrium (LD) analysis was conducted, and SNPs were greedily pruned until no pairs with r2 > 0.2 within a 500 kb sliding window with a step size of 50 markers (note that this value was updated from the originally preregistered 50 kb window upon a reviewer request at the revision stage of the manuscript). The list of SNPs selected after the process is reported in Extended Data Figure 1-3. Population structure and kinship were modeled as described in the statistical analyses below.
Neuroimaging data preprocessing and parcellation
With the goal of investigating specifically those connectome phenotypes that we previously established as being heritable (Jun et al., 2022), all neuroimaging procedures correspond to those that were performed in our prior work and resulted in the preregistered phenotype selection. The procedures are included here for completeness. Imaging data were acquired on a customized Siemens 3 T Skyra at Washington University in St. Louis using a multiband sequence. Each subject had four 15 min resting-state fMRI runs. We used the data that were minimally preprocessed by the HCP consortium (Glasser et al., 2013) with the pipeline of Smith et al. (2013a,b) using tools from FSL (Jenkinson et al., 2012) and FreeSurfer (Fischl, 2012) and artifact removal using ICA + FIX (Griffanti et al., 2014; Salimi-Khorshidi et al., 2014). Intersubject registration of the cerebral cortex was carried out using areal feature-based alignment and the Multimodal Surface Matching algorithm (“MSMAll”; Robinson et al., 2014; Glasser et al., 2016). The HCP team had previously parcellated the neuroimaging data using Independent Component Analysis (ICA) in FSL with model orders of 25, 50, 100, 200, and 300 and provided averaged BOLD time series data for regions of these group ICA-based parcellations. For the present study and our previous work (Jun et al., 2022), the 300-model order ICA was chosen, and its time series data was used without global signal regression. However, Automated Anatomical Labeling (AAL) atlas-based masks and visual inspection identified that 161 out of 300 independent components were either from cerebellar or brainstem. Therefore, the focus of our investigation was on the 139 independent components from the cortical and subcortical cerebrum, which represent the canonical intrinsic connectivity networks. The reconfigurations of these canonical networks have been described as key contributors to the behavioral impact of connectome dynamics (Thompson et al., 2013; Sadaghiani et al., 2015; Douw et al., 2016). Further details on cerebral parcellation can be found in our previous work (Jun et al., 2022).
Hidden Markov modeling and temporal features of the connectome dynamics
The following briefly describes the procedures performed in our prior twin heritability study that provided the foundation for the current investigation (Jun et al., 2022). The hidden Markov model (HMM) assumes that time series data can be represented by a finite sequence of hidden states. Each HMM-inferred connectome state, along with its corresponding time series, represents a unique connectivity pattern that temporally reoccurs over time. Using the HMM-MAR (multivariate autoregressive) toolbox (Vidaurre et al., 2016), we applied the HMM to the minimally preprocessed, region-wise BOLD time series and obtained four discrete connectome states (K of 4; Extended Data Fig. 1-4). Specifically, the number of states (K) was chosen following our previous work (Jun et al., 2022) which identified heritability of the temporal connectome dynamics phenotypes using K = 4 and replicated the findings using K = 6. Such Ks were selected based on prior HMM studies that have directly compared several Ks and identified Ks between 3 and 7 as optimal (Vidaurre et al., 2016; Stevner et al., 2019). Further details on the HMM and reasonings for the selection of K can be found in our previous work (Jun et al., 2022).
HMM-derived estimates describe the temporal aspects of connectome dynamics by characterizing the sequence of connectome states or the connectome's trajectory through state space. Specifically, we calculated Fractional Occupancy (the cumulative total time spent in a given state; 1 × K) and Transition Probability (the probability matrix of transitioning between all possible pairs of discrete states; K × K) for each subject. While Transition Probability and Fractional Occupancy are not fully independent measures, they do provide nonoverlapping information about connectome dynamics. For example, a state with particularly high Fractional Occupancy is likely to have high values as initial state and target state in the Transition Probability matrix. Despite such dependence, however, two hypothetical subjects with highly comparable Fractional Occupancy values across the k states may still have substantially different Transition Probability matrices. Our previous work demonstrated strong genetic effects on these temporal characteristics of functional connectome reconfigurations (Jun et al., 2022).
Our phenotypes of interest (i.e., Fractional Occupancies and Transition Probabilities) span multiple brain states and are inherently multidimensional. While recent methods extend genetic association analyses to multivariate phenotypes, these approaches involve several assumptions and may increase type I error rates (Chung et al., 2019). As our prior work (Jun et al., 2022) demonstrates the heritability of the multivariate phenotypes after dimensionality reduction, we believe that important phenotypic information is retained after such reduction.
To reduce the dimensionality of Fractional Occupancy and Transition Probabilities as in our prior work (Jun et al., 2022), we calculated the distance of individuals from an origin point. Specifically, we generated surrogate data to derive the origin point from which the Euclidean distance of each multivariate phenotype was estimated for each subject. Surrogate datasets were simulated using the simhmmmar function from the HMM-MAR toolbox, which preserves the static covariance structure of the original data while disrupting the precise temporal ordering of states. An HMM inference with K = 4 states was applied to each of 50 surrogate datasets, and Fractional Occupancy and Transition Probability were recomputed. The Fractional Occupancy (respectively Transition Probability) values averaged across the surrogate datasets were used to determine the origin point. Further details on obtaining the origin point from the null models can be in our previous work (Jun et al., 2022).
Cognitive performance
The following briefly describes the factorization that was performed in our prior heritability study and resulted in preregistered cognitive phenotypes of interest (Jun et al., 2022). Specifically, we included 14 cognitive measures provided by the HCP, which are summary scores from cognitive tasks or questionnaires under the cognition domain (see Extended Data Fig. 1-5 for more information on each variable). The measures were normalized to have zero mean and unit variance. A factor analysis was then conducted to cluster the cognitive measures into four Cognitive Factors (Jun et al., 2022). The Cognitive Factors were interpreted as follows: factor 1: “Language”; factor 2: “Impulsivity/self-regulation”; factor 3: “Cognitive control”; and factor 4: “Memory”. Further details on the factor analysis can be found in our previous work (Jun et al., 2022).
Statistical analysis
We evaluated the joint effect of the neuromodulatory SNPs on the temporal phenotypes of the connectome dynamics, i.e., Fractional Occupancy and Transition Probability. Specifically, separately for each neurotransmitter system, we aimed to build a model incorporating the optimal combination of SNPs and covariates that best explains each temporal connectome dynamics phenotype. Note that SNP genotypes were dummy coded as 0/1/2 to represent subjects with 0, 1, or 2 copies of the minor allele, which is according to the standard additive model applied in genetic association studies.
To avoid overfitting, we employed a fivefold cross-validation approach (Fig. 1C; Varoquaux et al., 2017). This approach involves splitting the full data into equally sized five folds and using four-fifths of the data as train set and the remaining one-fifth of the data as test set, then rotating which one-fifth is used for testing. When combined, the five test sets form the original full data. It is important to note that participants from the same family were assigned to the same fold.
In each train set, we used the Stability Selection procedure (Meinshausen and Bühlmann, 2010; Shah and Samworth, 2013) with elastic net regression as the selection algorithm to identify a subset of predictors that best explains each of the connectome phenotypes. Specifically, for each connectome phenotype, SNPs of a modulatory neurotransmitter system and covariates [i.e., age, sex, frame-wise displacement (FD) calculated across the entire scan duration for head motion, and 10 principal components that account for European-ancestry population structure] entered a Stability Selection procedure (Meinshausen and Bühlmann, 2010; Shah and Samworth, 2013) through the R package stabs (https://cran.r-project.org/web/packages/stabs/index.html). We opted for Stability Selection procedure as our selection algorithm because it offers greater stability compared with other regularization methods (e.g., stepwise regression or elastic net without Stability Selection; Hofner et al., 2015). The method also generates more parsimonious models by reducing the unmodeled variability to the estimates and, thus, mitigating the risk of type 1 errors. As such, Stability Selection has been widely used in the fields dealing with high-dimensional data, such as gene regulatory networks (Haury et al., 2012), genome-wide association studies (Patin et al., 2018), and graphical modeling (Bühlmann et al., 2014; Drton and Maathuis, 2017).
Specifically, Stability Selection implemented a complementary pairs subsampling approach (Shah and Samworth, 2013; Hofner et al., 2015), which involved dividing a given train set into 50 random pairs of subsamples (where each subsample comprised half of the data) and conducting variable selection on each subsample. Variables selected in a predetermined proportion of the subsamples were retained. Typically, the proportion threshold is preselected in the recommended range of 0.6–0.9 (Meinshausen and Bühlmann, 2010). We chose a threshold of 0.75 that balances rigorous selection with the risk of false negatives (Bommert et al., 2022).
Then, we fitted two models to the values of a connectome phenotype in the test set (i.e., the left-out fold): (1) the “full model” comprising a subset of covariates and SNPs selected by Stability Selection on the train set, along with random genetic effects in the form of a kinship matrix, and (2) the “base model,” which included only the selected covariates and the kinship matrix (excluding SNPs). Note that the base model was identical to the full model, except for the absence of SNPs. The kinship matrix, representing family structure, was obtained using the lmekin function in the R package coxme (https://cran.r-project.org/web/packages/coxme/index.html). Importantly, we did not predict the test set using the parameter weights of these two models derived from the training set. Instead, we fitted the two models (i.e., full and base models) directly to the test set and obtained the residuals, respectively. This approach was necessary because accounting for the kinship in the linear mixed-effects models, a critical step to adjust for the random genetics effect in the data with families, does not allow using the predefined parameter weights obtained from the training set and, obviously the kinship matrix differs between the training and testing samples. This process was repeated five times, with each fold serving as the test set once.
Consequently, each subject had two residuals: one from the full model and the other from the base model. The absolute values of these residuals were then concatenated across subjects, i.e., across the five test sets (Indahl and Naes, 1998; Eriksson et al., 2008). To assess whether the full model fits the data significantly better than the base model, the two sets of residuals entered an F test. In the regression context, the numerator degrees of freedom (df1) are related to the number of parameters being tested, which is the difference in the number of parameters between the full and the base models. The denominator degrees of freedom (df2) are related to the number of observations minus the number of parameters in the full model. In order to compute the F statistics and p values with the most conservative approach, we set the df1 to be largest by taking the maximum number of parameters from the full models and the minimum number of parameters from the base models across the five folds. All statistical analyses were performed in R.
Subsequently, in a second analysis, the above-described neurotransmitter models—specifically built for each of the connectome dynamics phenotypes—underwent further testing to evaluate their explanatory capacity for cognitive abilities. Specifically, each of the neurotransmitter models that outperformed the base model was tested for a new set of phenotypes: the Cognitive Factors, i.e., Language, Impulsivity, Cognitive Control, and Memory. The “cognitive full model” was fitted to the test set, which includes the same subset of SNPs and covariates selected above for each fold and the kinship matrix as predictors and each of the Cognitive Factors as a dependent variable. Simultaneously, the “cognitive base model” was fitted to the same test set, consisting only of the above-selected covariates and the kinship matrix. The significance of each cognitive model against their respective cognitive base model was tested using the same approach (Extended Data Table 2-1).
Table 2-1
Neurotransmitter models concomitantly predictive of Cognitive Factors. Modulatory neurotransmitter SNPs and covariates selected using stability selection algorithm on the train set (4/5th of the data) in each fold for the temporal connectome dynamics phenotypes were re-fitted to each of the Cognitive Factors. NT; Neurotransmitter, ACh; acetylcholinergic, NAd; noradrenergic, DA; dopaminergic, 5HT; serotonergic, PC; principal component that corrects for population structure, P†; P values Bonferroni-corrected for 8 tests. *P† < .05. Download Table 2-1, DOCX file.
Results
Consequent to the selection and quality control process of the genotype data (Fig. 1A), a total of 155 SNPs were included in the present study (Extended Data Fig. 1-3): 25 SNPs from 14 ACh genes, 9 SNPs from 4 NAd genes, 38 SNPs from 5 DA genes, and 83 SNPs from 17 5-HT genes. Two heritable connectome phenotypes and four associated cognitive phenotypes were adopted from our prior twin study (Jun et al., 2022).
Effects of select SNPs on temporal connectome phenotypes
Stability Selection (Fig. 1C) identified a subset of SNPs for each neurotransmitter system that best explained the connectome states’ Fractional Occupancy or Transition Probability, respectively. The SNPs and covariates included for each model by Stability Selection are listed in Extended Data Table 1-1. We found two neurotransmitter models that significantly outperformed the covariates-only base model (p† < 0.05, where p† is the Bonferroni-corrected p value for eight models; Table 1). Specifically, Fractional Occupancy was explained by the noradrenergic model (F(2,671) = 6.56, p† = 0.012), and Transition Probability was explained by the cholinergic model (F(3,670) = 12.25, p† = 6.61 × 10−7). The significant models included SNPs related to two NAd alpha 1 receptors, both G-protein-coupled receptors (ADRA1A and ADRA1D genes), an ACh metabolic enzyme (BCHE gene), and ACh nicotinic ion channels (CHRNB1, CHRNB4, and CHRNA2 genes). To provide an intuition, Figure 2 illustrates how one polymorphism contributing to a significant association model is linked to each connectome dynamics phenotype. Note that while the illustration highlights a single SNP for demonstration purposes, the significance of the models emerged from the synergistic impact of the respective SNPs in combination.
An illustration of genetic polymorphisms associated with temporal connectome dynamics phenotypes. To provide an intuition of the association, the figure uses the example of rs35105284, an SNP included in the significant noradrenergic (NAd) model of Fractional Occupancy. Averaged across the subjects of a given genotype, the pie charts depict the mean duration of time spent in each connectome state. For example, subjects with the homozygous genotype for major allele (T/T) of rs35105284 spend more time in state 4 and less time in state 3 compared with minor allele carriers (T/C and C/C). Note that while this illustration highlights a single SNP for demonstration purposes, the significance of the models often emerged from the synergistic impact of several SNPs in combination.
Neurotransmitter models explaining temporal connectome dynamics phenotypes
Table 1-1
Predictors selected from stability selection across five training sets for each temporal connectome dynamics phenotype and neurotransmitter system. Note that the variant type information of each variant (or SNP) has been acquired from (https://www.ensembl.org; accessed on 24 May 2023). ACh; acetylcholinergic, NAd; noradrenergic, DA; dopaminergic, 5HT; serotonergic, GPCR; G-protein coupled receptor, 5-HT; 5-hydroxytryptamine. Download Table 1-1, DOCX file.
Table 1-2
Gene expression levels of the select genes in the seven canonical intrinsic connectivity networks (ICNs). [A] A glass brain view of gene expression levels for ADRA1D and BCHE as examples of selected genes. [B] Bar plots depict the gene expression levels of all selected genes as listed in Table 1 (with exception of CHRNB4, which was not available in the microarray data used in this analysis). The gene expression level was averaged by the ICNs as defined by (Yeo et al., 2011). Expression levels were obtained from the Allen Institute for Brain Science (http://human.brain-map.org/static/download; accessed on Oct 12, 2022). For information regarding the normalization process of the microarray data and the units of intensity used in this dataset, see http://help.brain-map.org/download/attachments/2818165/Normalization_WhitePaper.pdf?version=1&modificationDate=1361836502191&api = v2. Download Table 1-2, TIF file.
For completeness, we also tested each of the 155 SNPs individually (along with the covariates) in separate linear regression models. This analysis represents our originally preregistered association approach prior to identifying and adopting an established and stable approach (Meinshausen and Bühlmann, 2010; Shah and Samworth, 2013) to unveil the impact of multiple SNPs. When tested individually in this manner, none of the SNPs survived a highly conservative correction method (Bonferroni-corrected p < 0.05/155), strengthening the notion that the SNPs are associated collectively, rather than individually, with the investigated connectome dynamics phenotypes.
Effects of select SNPs on cognition
In accord with our preregistered hypotheses, we tested whether the above-described models containing a set of SNPs that explain connectome dynamics phenotypes can also explain cognitive abilities. Specifically, the NAd model of Fractional Occupancy and the ACh model of Transition Probability were fitted to each of the four Cognitive Factors. In other words, the multi-SNP models initially built for connectome phenotypes were refitted to cognitive phenotypes. We found two models explaining cognitive abilities compared with their respective base model (p† < 0.05, where p† is the Bonferroni-corrected p value for 8 models in Table 2; see full results in Extended data Table 2-1). Specifically, the NAd SNPs originally identified in each fold to explain Fractional Occupancy explained Language and Memory factors.
Neurotransmitter models explaining cognitive factors
Discussion
The complexity of brain function and cognitive abilities arise from brain dynamics, including the functional connectome that changes dynamically over time. Such time-varying connectome dynamics, which are highly heritable, predictive of cognitive abilities (Vidaurre et al., 2017; Jun et al., 2022), and linked to mental health (Wang et al., 2020; Chen et al., 2023), may be modulated by the modulatory neurotransmitter systems (Shine et al., 2018b, 2019; Lurie et al., 2020). However, the specific genetic polymorphisms that shape connectome dynamics have been largely unknown. The current study departs from previous genetics neuroimaging investigations of static functional connectome measures (Whelan et al., 2012; Markett et al., 2016; Simola et al., 2022; Wang et al., 2022) and reveals that genetic profiles of modulatory neurotransmitter systems impact dynamic connectome phenotypes, specifically Fractional Occupancy and Transition Probabilities. Our findings provide evidence for genetic effects on connectome dynamics via the regulatory impact of modulatory neurotransmitter systems.
Interestingly, the genetic impact on dynamic connectome phenotypes was observed in the NAd and ACh systems, the two modulatory neurotransmitter systems that have been hypothesized to drive functional connectome dynamics through their cooperative and competitive impact on neural gain (output as a function of input; Shine, 2019). Notably, we observed a double dissociation in how these neurotransmitter systems contribute to the connectome phenotypes, with the NAd system (specifically alpha 1 receptor genes) explaining Fractional Occupancy and the ACh system (most notably nicotinic ion channel and enzyme genes) explaining Transition Probability. This observation indicates that distinct neural mechanisms support Transition Probability and Fractional Occupancy, which contain nonoverlapping information about connectome dynamics (although they are not fully independent measures; compare Materials and Methods). However, understanding of the mechanisms behind the differential effects of modulatory neurotransmitters on connectome phenotypes will require dedicated studies in the future.
The observed contribution of neurotransmitter systems to the heritability of connectome dynamics aligns with the role of these systems in driving widely distributed but spatially organized neuromodulation. Neuromodulatory neurotransmitter receptors and transporters are spatially distributed according to canonical intrinsic connectivity networks and functional gradients that are hallmarks of the connectome's architecture, thereby effectively mediating the connectome-wide propagation of synaptic communication and population activity (Shine, 2019; Goulas et al., 2021; Hansen et al., 2022). Extended Data Table 1-2 showcases this distributed gene expression across the intrinsic neurocognitive networks (Yeo et al., 2011) for each of the genes identified in Table 1 using the microarray data of the Allen Institute for Brain Science. In dedicated studies, this spatial organization may provide mechanistic insights into how genetic variation in neurotransmitter systems may modulate connectome states’ Fractional Occupancy and Transition Probabilities. As an example in prior literature, combining spatial mapping of a serotonergic receptor with modeling and pharmacological intervention (serotonergic agonist psilocybin), the causal chain from serotonergic neuromodulation to connectome state's Fractional Occupancy and Transition Probabilities has been outlined (Lord et al., 2019; Kringelbach et al., 2020). Building on our current findings, such models can be extended to include the subjects’ receptor/enzyme/transporter genotypes to mechanistically understand interindividual variability in connectome-level effects of neuromodulation.
The identified genetic impact of modulatory neurotransmitter systems on connectome dynamics may help understand the patho-etiology of various neurological and psychiatric conditions and predicting treatment response. For instance, the ADRA1D gene was consistently identified in the noradrenergic model of Fractional Occupancy (Table 1, Extended Data Table 1-1). ADRA1D, along with ADRA1A which contained all other identified NAd-related SNPs, encode NAd alpha 1 receptors whose activation contributes to cortical arousal (including enhanced excitability Datta et al., 2019; Porter-Stransky et al., 2019). Alpha 1 receptors thus serve as a pharmacological target (for receptors-specific antagonists) in posttraumatic stress disorder characterized by hyperarousal (Datta et al., 2019). The observed association between ADRA1 receptors and connectome dynamics may help explain the previously established link between physiological arousal fluctuations and functional connectome dynamics (Raut et al., 2021). Further, the BCHE gene was repeatedly selected in the cholinergic model of Transition Probability (Table 1, Extended Data Table 1-1). The BCHE gene encodes the butyrylcholinesterase enzyme that inactivates acetylcholine and consequently its modulatory impact on processing of environmental stimuli (Picciotto et al., 2012). BCHE is thought to contribute to accelerated deactivation of acetylcholine in the aging brain (Darvesh, 2016; Jasiecki et al., 2021) and is implicated in mild cognitive impairment (Li et al., 2021) and Alzheimer's disease (Darvesh, 2016; Jasiecki et al., 2021) While BCHE expression in astrocytes is widely distributed across the brain, the expression in neurons is particularly high in limbic areas (Darvesh et al., 1998; Darvesh and Hopkins, 2003; as confirmed in Extended Data Table 1-2), regions that play a critical role in Alzheimer's dementia. Butyrylcholinesterase is considered a potential therapeutic target for Alzheimer's disease, which is characterized by a cholinergic deficit (Diamant et al., 2006). It is interesting to note that the specific rs1803274 variant of BCHE gene found in our study has been strongly associated with the serum level of the butyrylcholinesterase (Benyamin et al., 2011). As connectome features have the potential to inform individualized medicine (Damoiseaux et al., 2021), the genes we identified may help to refine treatment efforts that embrace interindividual differences in connectome dynamics.
In addition, we found that the noradrenergic model, which explained Fractional Occupancy, also explains cognitive factors characterized by performance in the domains of language and memory (Table 2). It is important to note that the specific emergence of Language and Memory factors in this association closely aligns with our previous heritability work in the same dataset (Jun et al., 2022). Specifically, the prior work showed that connectome phenotypes that were identified as heritable (and hence entered the current study) exhibit the strongest cognitive association with the same Language factor, followed by the Memory factor. Interestingly, a prior study (Han and Adolphs, 2020) that investigated heritability of behavioral measures, likewise in the HCP dataset, demonstrated the highest heritability in two cognitive factors with task loadings similar to those of our Language and Memory factors. Collectively, these findings not only suggest that both subject-specific trajectories across connectome states and cognitive abilities are under genetic effects but also that the noradrenergic polymorphisms may drive the relationships between connectome phenotypes and cognition (Shine et al., 2016, 2018a, 2019; Pfeffer et al., 2018; Totah et al., 2019), specifically between Fractional Occupancy and Language/Memory.
Our study is subject to several limitations which must be considered when interpreting our findings. Firstly, while we identified polymorphisms associated with connectome dynamics, in some cases it is yet to be revealed whether they change the functional capability of the respective proteins. Specifically, beyond missense protein-coding SNP the identified models contained introns and synonymous SNPs (e.g., rs35105284; ADRA1D gene). Such inclusion of synonymous SNPs and introns is in line with the observation that complex traits are primarily associated with noncoding variants that likely impact gene regulation (Welter et al., 2014). Furthermore, our work should be interpreted at the level of identified genes rather than individual SNPs, as we have focused on investigating variants which have significant influence on gene expression (eQTL p < 0.05), which is their most likely mechanism of action. Finally, it is important to note that the genes identified in our study and their encoded products (such as receptors or metabolic enzymes) function as components of intricate system-level networks (Vidal et al., 2011) influenced by various biological and environmental factors, making mechanistic interpretation of our results challenging. For example, different receptors have distinct binding affinities, downstream signaling pathways, and cellular responses to varying levels of neuromodulatory neurotransmitters. Additionally, neural activity can also impact gene expression and neurotransmitter function (Flavell and Greenberg, 2008; West and Greenberg, 2011). Thus, the associations between specific genes and connectome phenotypes identified in our study should be interpreted cautiously and further investigated in the future in terms of the underlying molecular and cellular mechanisms. The discoveries from our study motivate investigations into such potentially indirect (i.e., regulatory) impact of the identified SNPs on the respective proteins (Tam et al., 2019) and suggest that such impact translates into connectome states’ Fractional Occupancy and Transition Probability phenotypes. Moreover, our study is based on a relatively small size of the European-ancestry population due to the difficulty in finding a consortium-level dataset with publicly available genotype and high-quality resting-state fMRI data of young adults, calling for future replication upon availability of such data. Further, because ethnic differences in the frequency or the effect size of a variant can affect the likelihood of its discovery and its contribution to phenotypes across populations, a study cannot include subjects from multiple ethnic groups. Future studies should replicate our approaches in other ethnic groups and compare the roles of potentially different neurotransmitter genetic profiles in modulating the connectome dynamics and cognition.
To conclude, our study reveals specific genes that modulate dynamic connectome reconfigurations. This work extends beyond prior studies that have primarily focused on static connectome measures. We demonstrate that multiple SNPs in each neurotransmitter system are critical in shaping the dynamic trajectory of the functional connectome in state space. Our study has the potential to inform the development of precision medicine strategies that leverage genetic information to tailor diagnostics and interventions for mental health disorders.
Footnotes
This work was partly supported by the National Institute for Mental Health (1R01MH116226 to S.S.). Neuroimaging and behavioral data were provided by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University. The genotype data presented in the current publication are based on the use of study data downloaded from the dbGaP web site, under dbGaP Study Accession: phs001364.v1.p1.
The authors declare no competing financial interests.
- Correspondence should be addressed to Suhnyoung Jun at suhnyoung.jun{at}gmail.com.
