A developmental reduction of the excitation:inhibition ratio in association cortex during adolescence.

Adolescence is hypothesized to be a critical period for the development of association cortex. A reduction of the excitation:inhibition (E:I) ratio is a hallmark of critical period development; however, it has been unclear how to assess the development of the E:I ratio using noninvasive neuroimaging techniques. Here, we used pharmacological fMRI with a GABAergic benzodiazepine challenge to empirically generate a model of E:I ratio based on multivariate patterns of functional connectivity. In an independent sample of 879 youth (ages 8 to 22 years), this model predicted reductions in the E:I ratio during adolescence, which were specific to association cortex and related to psychopathology. These findings support hypothesized shifts in E:I balance of association cortices during a neurodevelopmental critical period in adolescence.

INTRODUCTION

Adolescent brain development is characterized, in part, by the continued structural and functional maturation of the association cortices (–). The specificity of the developmental timing and localization of association cortex maturation as well as the links between association cortex development and long-term psychiatric outcomes have led to the hypothesis that adolescence functions as a critical period of development within association cortex (). Critical periods are windows of development during which experience powerfully shapes the development of neural circuits through heightened experience-dependent plasticity with long-term impacts on behavior (). These important neurodevelopmental windows are theorized to progress hierarchically throughout development, beginning in primary sensory cortices and sequentially advancing to secondary and higher-order cortical areas (, ). The neurobiological mechanisms that underlie critical periods are thought to be conserved across the cortex and have been carefully delineated in decades of work on early critical periods in sensory cortex (–).

One of the hallmark features of critical period development is the maturation of GABAergic (γ-aminobutyric acid–mediated) inhibitory circuitry, particularly parvalbumin (PV)–positive interneurons, leading to a reduction in the excitation:inhibition (E:I) ratio (). The reduction of the E:I ratio leads to an increase in the signal-to-noise ratio of local circuit activity as inhibition suppresses the effect of spontaneous activity on circuit responses in favor of stimulus-evoked activity (). This essential mechanism has been shown to regulate the timing of critical period development across visual (), auditory (), and sensorimotor cortices (). As such, if the adolescent critical period hypothesis is correct, inhibitory maturation should result in a developmental reduction in the E:I ratio across adolescence within association cortex.

Evidence for E:I development in association cortex during adolescence has been largely limited to animal models. This work has suggested that prefrontal GABAergic inhibitory circuitry undergoes substantial modifications. Specifically, PV interneurons, a critical component of E:I maturation in sensory system critical periods, have been shown to increase in prefrontal cortex during adolescence in the rat () and nonhuman primate (). At the same time, the expression of GABAA receptor α1 subunits, which are primarily expressed on PV cells and support fast synaptic inhibition () and synaptic plasticity (), also increase during adolescence in the prefrontal cortex of the nonhuman primate (, ). These neurobiological changes lead to important functional increases in inhibitory signaling, effectively reducing the E:I ratio (, ). Together, these findings are suggestive of critical period development and may indicate that similar processes are unfolding in the human (). Translating these findings to human studies of development is crucial, as disruptions to the E:I balance are hypothesized to play a significant role in the onset of psychiatric disorders (, ). However, the extent to which these critical period mechanisms are present in association cortex during adolescence in the human remains largely unexplored. Corroborating evidence has been found in postmortem studies, which demonstrate increases in PV () and GABAA α1 expression (), but it has been unclear how to measure developmental changes in the E:I ratio in vivo in humans using available neuroimaging techniques. This lack of in vivo measures has limited our ability to test the adolescent critical period hypothesis.

Here, we leveraged a pharmacological functional magnetic resonance imaging (fMRI) (phMRI) experiment using a GABAergic benzodiazepine challenge to empirically generate a model for the effect of inhibitory modulation of the E:I ratio on patterns of fMRI connectivity. Benzodiazepines are positive allosteric modulators of the GABAA receptor that increase the effectiveness of postsynaptic GABAergic signaling, resulting in an increase in inhibition relative to excitation. Benzodiazepines have been used to pharmacologically manipulate the E:I ratio by enhancing inhibitory signaling in disease models of E:I imbalance (, ) as well as in studies of critical period development (, , ). In the current study, we first trained a multivariate model to distinguish benzodiazepine-induced change in the E:I ratio and established the neurobiological relevance of our empirical model by comparing the model features to known aspects of benzodiazepine pharmacology as well as a functional gradient that has been shown to reflect patterns of excitatory and inhibitory interneuron expression (, ). We then applied our trained and validated model to a large independent developmental dataset to investigate E:I changes occurring in association cortex during adolescence. We hypothesized that patterns of functional connectivity would develop to reflect a reduction in the E:I ratio that is specific to association cortex.

RESULTS

An empirical model of the E:I ratio

Forty-three adult participants completed a double-blind, placebo-controlled phMRI study with the benzodiazepine alprazolam (86 sessions total). Alprazolam is a classical benzodiazepine that enhances the effect of GABA at GABAA receptors through positive allosteric modulation, increasing inhibition and effectively reducing the E:I ratio (). Participants were orally administered 1 mg of alprazolam, which produces an increase in GABAergic inhibition that is considered to be clinically effective (). Functional connectivity matrices were derived for placebo and drug phMRI sessions using a top-performing pipeline that minimized the impact of motion artifact (). A linear support vector machine (SVM) classifier was trained to distinguish placebo and drug sessions based on the multivariate patterns of functional connectivity (Fig. 1, green pathway). Cross-validation and permutation testing revealed that the trained SVM identified drug versus placebo sessions in left-out data far better than chance [area under the receiver operating curve (AUC) = 0.716, ppermutation = 0.002; Fig. 2A]. Classification accuracy was not significantly associated with participant age, participant sex, drug condition (i.e., accuracy did not differ between drug and placebo sessions), session order (day 1 versus day 2), self-reported feeling of relaxation, blood concentration of alprazolam following the experiment, state anxiety before or after scan, trait anxiety, or head motion during the scan (table S1).

Fig. 1.

Analysis workflow.

Dataset: Two datasets were collected on the same scanner using highly similar acquisition parameters: a phMRI dataset using the benzodiazepine alprazolam (green) and a developmental fMRI sample from the Philadelphia Neurodevelopmental Cohort (PNC) (purple). Preprocessing: Datasets were preprocessed using identical pipelines, which included removal of nuisance signal with aCompCor, global signal regression, and task regression. Connectivity matrix generation: Connectivity matrices were generated from standard atlases for placebo and drug sessions from the alprazolam dataset (n = 43; 86 sessions total) and for the PNC dataset (n = 879). Train and validate model: The alprazolam dataset was used to train a linear SVM classifier to distinguish drug and placebo sessions using 10-fold cross-validation. Apply model: The validated alprazolam model was applied to the PNC dataset, generating a distance metric that reflected each participant’s position on a continuum from “drug-like” (lower E:I) to “placebo-like” (higher E:I). Regress model output on age: This metric was then regressed on age using a generalized additive model with penalized splines that included covariates for sex, head motion, and attentiveness.

Fig. 2.

A multivariate model distinguishes alprazolam and placebo sessions, capturing E:I ratio.

(A) Classifier performance. The binary SVM classifier identified drug and placebo sessions in 10-fold cross-validation with an AUC of 0.716 and an accuracy of 69.5% (top). The observed AUC and accuracy were significantly greater than a permuted null distribution (bottom). (B) Mean absolute feature weights for all nodes from the validated SVM model.

Furthermore, sensitivity analyses confirmed that in-scanner head motion was not associated with our pharmacological manipulation or model performance (fig. S1). The spatial pattern of estimated feature weights from the SVM model highlighted the contributions of subcortical regions, including the thalamus and amygdala, and also contributions throughout the cortex (Fig. 2B).

Biological relevance of the E:I model

Next, we established the biological relevance of the trained E:I model. First, we compared the spatial pattern of cortical feature weights to a widely used functional gradient of macroscale cortical organization that places regions on a continuum from unimodal to transmodal function (). This continuum has been shown to capture variation in excitatory neuron structure, inhibitory interneuron expression, and E:I balance (, ). Using a recently developed analytic procedure that accounts for spatial autocorrelation structure (), we observed a significant relationship between our model feature weights and this pattern of macroscale cortical organization (r = 0.33, P = 0.003; Fig. 3A). This finding suggests that GABAergic modulation of functional connectivity patterns varies along a transmodal-to-unimodal gradient that, in part, indexes diversity in excitatory and inhibitory neurobiological properties. Next, we investigated whether the estimated model features corresponded to the known pharmacology of benzodiazepines like alprazolam. Of the six GABAA subunit receptors, GABAA α1 to α6, only α1, α2, α3, and α5 are sensitive to benzodiazepines because of the presence of an amino acid residue, histidine (). We used the Allen Human Brain Atlas () to evaluate how the feature weights from the classifier model aligned with spatial patterns of gene expression for the six GABAA subunit receptors, GABRA1 to GABRA6 (corresponding to GABAA α1 to α6). We found evidence of a clear biological double dissociation: Model features were significantly associated with the gene expression patterns of the benzodiazepine-sensitive GABAA subunits (α1, α2, α3, and α5) and not the benzodiazepine-insensitive GABAA subunits (α4 and α6; Fig. 3B).

Fig. 3.

Model features align with cortical organization and benzodiazepine pharmacology.

(A) The cortical pattern of nodal SVM weights from the multivariate E:I ratio model was significantly associated with transmodality using an established measure of macroscale cortical organization (). (B) Nodal weights were also specifically correlated with the spatial patterns of benzodiazepine (BZD)–sensitive GABAA receptor subunit expression. Spatial relationships were tested for significance against a spatial autocorrelation–preserving null distribution [BrainSMASH ()] and corrected for multiple comparisons using the Bonferroni correction (PBonf). n.s., not significant.

Development of the E:I ratio during adolescence

We next used our empirically generated E:I ratio model to test the hypothesis that the E:I ratio declines as part of the critical period of association cortex development. An independent sample of 879 youth (aged 8 to 21.7 years) participated in a highly similar fMRI acquisition on the same scanner; these data were preprocessed using an identical pipeline. We applied our validated E:I model to the developmental dataset without further tuning and obtained the model-estimated distance from the classification hyperplane. This metric reflects a participant’s position on the continuum between “drug-like” (lower E:I) and “placebo-like” (higher E:I). To capture both linear and nonlinear effects in a rigorous statistical framework, we then regressed this metric on age using a generalized additive model (GAM) with penalized splines (Fig. 1, purple pathway). We found that age was positively associated with patterns of GABA-modulated functional connectivity, reflecting an age-related reduction in the E:I ratio. Significant reductions occurred between ages 12.9 and 16.7 years [Fs (Age) = 3.11, P = 0.037], after which the slope of the effect no longer significantly differed from zero, indicating a period of developmental stability (table S2). The age-related reduction in the E:I ratio was robust across multiple alternative parcellation schemes (table S3 and fig. S2).

Age-related reductions in the E:I ratio are specific to association cortex

We hypothesized that age-related reductions in the E:I ratio during adolescence were specific to association cortices. To test this hypothesis, we trained two additional models that restricted input features to connections to the most transmodal (Fig. 4A, blue) or unimodal (Fig. 4A, green) parts of the cortex. Both models significantly distinguished drug from placebo phMRI sessions (Fig. 4A). However, when applied to the developmental dataset, significant age-related reductions in the E:I ratio were only observed for the model trained on connections with transmodal cortex [transmodal: Fs (Age) = 9.96, P = 0.0017; unimodal: Fs (Age) = 3.59, P = 0.058; transmodal versus unimodal: Fs (Age) = 5.96, P = 0.015; Fig. 4B]. These results suggest that transmodal association cortices undergo a reduction in the E:I ratio during adolescence, consistent with a critical period of development.

Fig. 4.

Transmodal areas undergo E:I ratio development during adolescence.

(A) Model performance for unimodal and transmodal classifiers. SVM classifiers were trained and validated for connections to the most transmodal (green) and most unimodal (blue) areas only. Dashed lines indicate acquisition field of view for the phMRI dataset. Both models performed significantly better than a permuted null distribution (middle: receiver operating characteristic curves for each model; right: null distributions from 1000 null permutations). (B) Models trained on transmodal and unimodal data were applied to the developmental dataset, generating a distance metric for each participant where greater values represent patterns of functional connectivity consistent with a lower E:I ratio. Individuals had a lower estimated E:I ratio with age in transmodal cortex (left) but not in unimodal cortex (middle). This pattern was confirmed by a significant effect of age on within-subject change in transmodal versus unimodal distance scores (right). *P < 0.05 and **P < 0.01.

Analysis of dimensions of psychopathology

Last, we investigated whether individual differences in dimensions of psychopathology (, ) were associated with the E:I ratio of association cortex. We found that mood disorder symptomatology, but not other psychopathology dimensions, moderated age-related differences in the estimated transmodal E:I ratio. Specifically, individuals with greater lifetime mood disorder symptoms displayed a relatively stable E:I ratio over development instead of the normative reduction of the E:I ratio (Age*Mood interaction: F = 7.64, P = 0.0058).

DISCUSSION

We used a unique combination of human phMRI and developmental fMRI data to provide evidence for an essential component of critical period development: developmental reductions in the E:I ratio. Our approach generated an empirical fMRI model of the E:I ratio that showed a high degree of correspondence to known GABAergic benzodiazepine neuropharmacology and that could be applied to a large independent sample of youth. Consistent with our hypothesis, this approach revealed that patterns of functional neurodevelopment in adolescence are consistent with developmental reductions in the E:I ratio that are specific to association cortex. Furthermore, we show that individual differences in this process are associated with individual differences in lifetime mood symptom burden, in alignment with models positing that E:I abnormalities underlie the emergence of psychopathology (, , , ). Together, these findings support the hypothesis that critical period mechanisms shape association cortices during adolescence.

Critical period development has been predominantly associated with early sensory cortex development. Since the first studies of critical period development in the visual cortex almost 60 years ago (), a wealth of previous work has elucidated the mechanisms that shape critical period plasticity in these areas. These studies have identified the maturation of local inhibitory circuitry, particularly PV interneurons, and its resulting impact on the E:I balance as an essential critical period mechanism (, ). This phenomenon is necessary for the opening of the critical period window, facilitates critical period plasticity, and is present in critical periods across sensory modalities (, , , ). For example, modulation of the E:I balance using benzodiazepines has been shown to be sufficient to accelerate the timing of critical period plasticity (, ). The results of this study suggest that this phenomenon also occurs in association cortex during human adolescence.

Our findings align with a growing literature characterizing inhibitory maturation during this developmental stage. Animal models and postmortem human studies have shown maturation of inhibitory neurobiology in the prefrontal cortex during adolescence, including increasing expression of PV interneurons and GABAA α1 receptor subunits (, , , ). These processes increase functional inhibition, reducing the E:I ratio and increasing the signal-to-noise ratio of circuit activity (, , ). Computational simulations have suggested that these maturations also facilitate high-frequency oscillatory capability (). This is consistent with human electroencephalography studies showing increased γ-band oscillatory power during adolescence (). Last, two recent magnetic resonance spectroscopy studies have shown increases in GABA levels relative to glutamate levels in frontal cortex during adolescence (, ). Although it is not possible to examine the functional effect of these changes on the E:I ratio using spectroscopy, these findings align with a model of developmental reduction in the E:I ratio during adolescence. This body of previous work coheres with the findings presented here and is consistent with a critical period model of adolescent association cortex development. Just as sensory critical period plasticity refines neural circuits underlying sensory processing, the critical period for association cortex may facilitate the plasticity of circuits that underlie the higher-order cognitive processes refined during adolescence and are thought to be dependent on association cortex (, ).

It should be noted that there are two classes of critical period mechanisms: facilitating factors, which open the critical period window and facilitate plasticity, and braking factors, which stabilize neural circuits and physically limit future plasticity (, ). The maturation of inhibition and the resulting reduction in the E:I ratio are critical period facilitators (). Using GAMs, which can flexibly capture linear and nonlinear effects while penalizing overfitting, we found that the model fit for age-related reductions in the association cortex E:I ratio was linear. It is important to note that this does not necessarily mean that critical period plasticity is linearly increasing or that the critical period window is persistently open over the entire age range reported here. The developmental reduction in the E:I ratio is indicative of critical period opening, but it does not provide information about critical period closure. Closure of the adolescent critical period would be dependent on the development of braking factors, such as myelination and the formation of perineuronal nets (PNNs) (, ), which may follow distinct developmental trajectories. Consistent with a critical period model, many studies have provided evidence of myelination of association cortex and large white matter (WM) pathways linking association cortex to other areas of the brain that continues into adulthood, including histological (), myelin mapping (, ), and diffusion imaging (, ). At present, studies of PNN formation are limited to postmortem methods and animal models, which have demonstrated developmental increases in PNN formation in the prefrontal cortex from adolescence to adulthood (, ). Together, these studies indicate that critical period braking factors are forming during the transition from adolescence to adulthood, stabilizing neural circuits and closing the critical period window. However, to precisely demarcate the opening and closing of critical period plasticity during adolescence, future work that jointly investigates the developmental time course of critical period facilitators, such as the E:I ratio reported here, and critical period braking factors is needed.

Mood-related psychopathology typically first emerges during adolescence, with adolescent onset predicting greater illness duration and comorbidity (). Here, we observed that beginning in adolescence, youth with greater burden of mood symptoms exhibit an altered pattern of E:I development within the association cortex. Previous cross-species research has linked E:I disruptions to the manifestation of mood symptomatology, with alterations in both GABA and glutamate neurotransmitter systems being implicated (, ). Animal models of depression exhibit lower expression of inhibitory interneuron markers and GABA synthesis enzymes, as well as reduced excitatory neuron dendritic spine and synapse number (). Human studies have provided convergent evidence, demonstrating reduced GABA levels in the brain in those with depression (). Our study supports E:I-based models of mood psychopathology and further places them within a key neurodevelopmental framework—underscoring how E:I disruptions may reflect atypical critical period development. As alprazolam modulates the E:I ratio via GABA, our findings suggest that mood symptoms in youth may, in part, be linked to impaired development of inhibition in transmodal cortex that supports social cognition, emotion regulation, and reward processing ().

Last, we note that the approach used in this study highlights the potential for phMRI data to generate insights into independent datasets to inform new hypotheses. We combined machine learning and phMRI using a GABAergic alprazolam challenge to generate an empirical model for the effect of GABAergic modulation on patterns of fMRI connectivity. As evidence for the efficacy of this approach, the trained model could not only significantly predict drug versus placebo sessions in unseen data but also the model features demonstrated a significant correspondence with known benzodiazepine neuropharmacology. Notably, the model features were significantly associated with the GABA receptor most strongly implicated in critical period development, the GABAA α1 receptor. Although, because of the rarity of phMRI data, we were not able to confirm the generalizability of our model with an independent alprazolam phMRI dataset, the model performance and underlying interpretability of the learned features highlight the biological relevance of this method. Whereas in this study we applied this method to an independent developmental dataset to provide insights into the critical period mechanisms unfolding during adolescence, future work could apply this approach to other datasets to inform new research questions.

Despite the strengths of this study, two limitations should be noted. First, recent work has demonstrated that functional connectivity estimates become more stable and reliable as the duration of the time series increases, suggesting that time series greater than 30 min may be ideal (, ). The datasets used in the current study acquired data that were considerably shorter in duration. As such, future work could strengthen the results presented here by collecting longer functional acquisitions that would improve the stability of functional connectivity estimation and potentially improve the performance of the classifier. Second, the existing dataset we used to generate our empirical E:I model consisted of an adult sample that was considerably older than the developmental sample to which we applied the validated model. An advantage of our approach is that our model was trained in a sample in which inhibitory maturation has completed, allowing us to isolate the effect of the GABAergic modulation of functional connectivity patterns irrespective of any developmental effects. Although it may also be informative to train the model in a sample of comparable age, it would be ethically and practically difficult to conduct a pharmacological imaging study in a pediatric population. As such, the use of an adult sample is likely the best feasible approach.

Together, the findings of this study support the hypothesis that critical period mechanisms, such as the inhibition-induced reduction of the E:I ratio, shape association cortices during adolescence. Studying development from a critical period perspective provides a powerful mechanistic framework for understanding how experience and neurobiology interact to shape long-term cognitive, social, and psychiatric outcomes. A critical period model of adolescent development can draw on the history of detailed work on sensory critical periods to generate testable hypotheses for the mechanisms unfolding during adolescence in association cortex. Understanding these mechanisms are a necessary prerequisite to understanding of how experience, environment, and neurobiology contribute to differing neurodevelopmental trajectories in health and mental illness. This work thus lays the groundwork for future studies of the unique impact of experience on neurodevelopment and also suggests the possibility of targeted interventions during this critical window of vulnerability to psychopathology ().

MATERIALS AND METHODS

Participants and experimental procedures

Alprazolam sample

The alprazolam sample and study procedures have been described in detail in our earlier work (). Briefly, 47 adults participated in a double-blind, placebo-controlled pharmacological imaging study using the benzodiazepine alprazolam. Exclusion criteria have been previously described (), but we note that participants were excluded for current or recent benzodiazepine use, current or recent treatment with drugs known to affect benzodiazepine action, history of alcohol abuse or dependence, and history of psychiatric or neurological disorders. Each participant completed two identical experimental sessions approximately 1 week apart. In one session, participants were given a 1-mg dose of alprazolam, and in the other, they were given an identical appearing placebo. One milligram of alprazolam produces an increase in GABAergic inhibition that is considered to be clinically effective (). The order of administration was counterbalanced across participants. Alprazolam or placebo was administered 1 hour before the fMRI acquisition so that alprazolam levels and effects were near their peak at the time of data collection (). Following the MRI sessions, blood was drawn to measure alprazolam plasma levels (high-performance liquid chromatography assay with a lower limit of quantitation of 5 ng/ml performed by NMS Labs, Willow Grove, PA). During both sessions, participants completed an emotion identification task that lasted 10.5 min, while fMRI data were collected. Task-related fMRI results have been previously reported (). Four participants were excluded because of excess head motion in at least one session (see below) for a final sample of 43 participants and 86 sessions total (ages 20.9 to 59.4; M = 40.3, SD = 13.12, male/female = 24/19). Study procedures were approved by the University of Pennsylvania Institutional Review Board, and all participants provided written informed consent.

Developmental sample

Neuroimaging data were obtained from a community-based sample of 1476 youth (ages 8 to 21.9; M = 14.63; SD = 3.43; male/female = 698/778) that were part of the Philadelphia Neurodevelopmental Cohort (PNC). Data collection procedures and sample characteristics have been previously described in detail (, ). fMRI data were collected while participants performed the same emotion identification task as the alprazolam sample; this is also described in previous work (). From this original sample, 306 participants were excluded on the basis of health criteria, including psychoactive medication use at the time of study, medical problems that could affect brain function, a history of psychiatric hospitalization, and gross structural brain abnormalities. A total of 234 participants were excluded from further analysis due to head motion (see below), and 56 were excluded for poor structural image quality. In sum, following health exclusions and rigorous quality assurance, we retained 879 participants (ages 8.0 to 21.7 at first visit; M = 14.95; SD = 3.24; male/female = 383/496).

Neuroimaging acquisition

All data were collected on a Siemens Trio 3T as previously reported (). Whole-brain structural data were obtained with a 5-min magnetization-prepared, rapid acquisition gradient-echo T1-weighted (T1w) image (MPRAGE) using the following parameters: repetition time (TR) = 1620 ms; echo time (TE) = 3.87 ms; field of view (FOV) = 180 × 240 mm; matrix = 192 × 256; effective voxel resolution = 1 × 1 × 1 mm3. Blood oxygen level–dependent (BOLD) fMRI data were obtained as a slab single-shot gradient-echo echoplanar imaging (EPI) sequence using the following parameters: TR = 3000; TE = 32 ms; flip angle = 90°; FOV = 240 mm; matrix = 128 × 128; slice thickness/gap = 2/0 mm; 30 slices; effective voxel resolution = 1.875 × 1.875 × 2 mm3; 210 volumes. As previously described (), data were acquired in a FOV that included temporal, inferior frontal, and visual cortices as well as subcortical structures (Fig. 3A, gray boxes).

All neuroimaging data were collected on the same Siemens Trio 3T scanner as was used for the alprazolam dataset. The neuroimaging procedures and acquisition parameters have been previously described in detail (). Briefly, structural MRI was acquired with a 5-min MPRAGE T1w image [TR = 1810 ms; TE = 3.51 ms; inversion time (TI) = 1100 ms; FOV = 180 × 240 mm2; matrix = 192 × 256; effective voxel resolution = 0.9 × 0.9 × 1 mm3]. BOLD fMRI was acquired using similar acquisition parameters to the alprazolam dataset. BOLD fMRI scans were acquired as single-shot, interleaved multislice, GE-EPI sequence sensitive to BOLD contrast with the following parameters: TR = 3000 ms; TE = 32 ms; flip angle = 90°; FOV = 192 × 192 mm2 (whole-brain acquisition); matrix = 64 × 64; 46 slices; slice thickness/gap = 3/0 mm; effective voxel resolution = 3.0 × 3.0 × 3.0 mm3; 210 volumes.

Preprocessing of neuroimaging data

All preprocessing was performed using fMRIPrep 20.0.7 (RRID:SCR_016216) (), which is based on Nipype 1.4.2 () and XCP Engine (PennBBL/xcpEngine: atlas in MNI2009 version 1.2.3; Zenodo: http://doi.org/10.5281/zenodo.4010846). The neuroimaging data from the alprazolam and developmental datasets were processed using identical pipelines as described below.

Anatomical data preprocessing

The T1w image was corrected for intensity nonuniformity with N4BiasFieldCorrection, distributed with ANTs 2.2.0 (), and used as T1w reference throughout the workflow. The T1w reference was then skull-stripped with a Nipype implementation of the antsBrainExtraction.sh workflow (from ANTs) using OASIS30ANTs as target template. Brain tissue segmentation of cerebrospinal fluid (CSF), WM, and gray matter was performed on the brain-extracted T1w using FAST in FSL 5.0.9 (). Volume-based spatial normalization to MNI2009c standard space was performed through nonlinear registration with antsRegistration (ANTs 2.2.0) using brain-extracted versions of both the T1w reference and the T1w template.

Functional data preprocessing

The alprazolam dataset consisted of two BOLD acquisitions per participant (drug and placebo session), which were preprocessed individually. The developmental dataset consisted of one BOLD acquisition per participant. All BOLD acquisitions were processed with the following steps. BOLD runs were first slice time–corrected using 3dTshift from AFNI 20160207 () and then motion-corrected using mcflirt (FSL 5.0.9). A field map was estimated on the basis of a phase difference map calculated with a dual-echo Gradient-recalled echo sequence, processed with a custom workflow of SDCFlows inspired by the epidewarp.fsl script and further improvements in Human Connectome Project (HCP) Pipelines (). The field map was then co-registered to the target EPI reference run and converted to a displacement field map with FSL’s fugue and other SDCflows tools. On the basis of the estimated susceptibility distortion, a corrected BOLD reference was calculated for a more accurate co-registration with the anatomical reference. The BOLD reference was then co-registered to the T1w reference using bbregister (FreeSurfer), which implements boundary-based registration (). Co-registration was configured with nine degrees of freedom to account for distortions remaining in the BOLD reference. Six head motion parameters (corresponding rotation and translation parameters) were estimated before any spatiotemporal filtering using mcflirt. Last, the motion-correcting transformations, field distortion correcting warp, BOLD-to-T1w transformation, and T1w-to-template (MNI) warp were concatenated and applied to the BOLD time series in a single step using antsApplyTransforms (ANTs) with Lanczos interpolation.

Confounding time series were calculated on the basis of the preprocessed BOLD data. The global signal was extracted within the whole-brain mask. In addition, a set of physiological regressors was extracted to allow for component-based noise correction (CompCor) (). Anatomical CompCor (aCompCor) principal components were estimated after high-pass filtering the preprocessed BOLD time series (using a discrete cosine filter with 128-s cutoff). The aCompCor components were calculated within the intersection of the aforementioned mask and the union of CSF and WM masks calculated in T1w space after their projection to the native space of each functional run (using the inverse BOLD-to-T1w transformation). Components were also calculated separately within the WM and CSF masks. In this study, for each aCompCor decomposition, the k components with the largest singular values were retained such that the retained components’ time series were sufficient to explain 50% of variance across the nuisance mask (CSF and WM). The remaining components were dropped from consideration. The head motion estimates calculated in the correction step were also placed within the corresponding confounds file. The confound time series derived from head motion estimates and global signals were expanded with the inclusion of temporal derivatives and quadratic terms for each ().

Subject-level time series analysis was carried out in XCP Engine using FILM (FMRIB’s Improved General Linear Model). All event conditions from the emotion identification task were modeled in the GLM as 5.5-s boxcars convolved with a canonical hemodynamic response function. Each of the five emotions (fear, sad, angry, happy, and neutral) was modeled as a separate regressor. The temporal derivatives and quadratic terms for each task condition as well as the confounding aCompCor, global signal, and motion time series described above were included as nuisance regressors. The nuisance regression pipeline used here has been shown to be a top-performing procedure for mitigating motion artifacts (). Consistent with our previous work, participants in the alprazolam dataset were excluded from future analyses if mean framewise displacement exceeded 0.5 mm in either session. A more stringent threshold of 0.3 mm was applied to the developmental dataset; head motion was also included as a covariate in all developmental models (see below).

Connectivity matrix generation

Fully preprocessed fMRI data were used to generate mean time series within a set of atlas-defined brain regions for each participant. Cortical regions were defined according to the Schaefer 400 parcel cortical atlas (). To accommodate the restricted FOV of the alprazolam BOLD acquisition, the atlas was masked such that only parcels with greater than 95% coverage were included in connectivity analyses. Subcortical regions were defined using the Automated Anatomical Labeling (AAL) atlas (). Subcortical areas included the left and right caudate, putamen, accumbens, pallidum, thalamus, amygdala, hippocampus, and parahippocampal area. These cortical and subcortical atlases were combined and used to generate mean time series for each region in each dataset. Functional connectivity was calculated as the correlation coefficient of the time series for each pair of regions (20,503 unique pairs). As part of sensitivity analyses, we repeated this process after defining cortical areas using Schaefer 200 parcellation (), the Multi-modal Parcellation atlas (), the Gordon cortical atlas (), or the AAL ().

Pharmacological classification analysis

We used a linear SVM to classify drug versus placebo sessions in the alprazolam dataset based on multivariate patterns of functional connectivity. Linear SVMs find a hyperplane to separate two classes of data by maximizing the margin between the closest points (the support vectors). SVMs were implemented in R () using the e1071 library and were trained using a linear kernel and the default parameters. Model performance was evaluated using 10-fold cross-validation, iteratively selecting data from 90% of participants as training data and testing the trained model on data from the remaining 10% of participants. Across testing sets, the prediction accuracy and AUC were calculated to evaluate model performance. To ensure that our results were not driven by a specific cross-validation split, we repeated the entire 10-fold cross-validation procedure 100 times, drawing the 10-fold subsets at random each time. Performance metrics were finally averaged across the 100 iterations of the cross-validation procedure.

To evaluate whether model performance (i.e., the accuracy and the AUC) was significantly better than expected by chance, we performed a permutation test. Specifically, we reapplied the cross-validation procedure 1000 times, each time permuting the session labels (drug and placebo) across the training samples without replacement. Significance was determined by ranking the actual prediction accuracy versus the permuted distribution; the P value of the accuracy and AUC was calculated as the proportion of permutations that showed a higher value than the observed value in the real, unpermuted data. To confirm that correct versus incorrect classification of drug and placebo sessions were not influenced by demographic or experimental factors, we conducted additional analyses. Demographic variables included participant age, participant sex, total score on the Structured Interview for Schizotypy (SIS), and state (before and after scan) and trait anxiety scores from the State Trait Anxiety Inventory (STAI). Experimental factors included drug condition (whether accuracy differed between alprazolam versus placebo sessions), order of testing (whether accuracy differed between sessions on day 1 versus day 2), the blood plasma concentration of alprazolam following the experiment, self-reported feelings of relaxation during the experiment, and head motion during the fMRI scan (mean framewise displacement). Continuous variables were evaluated using independent samples t tests; categorical variables were evaluated using χ2 tests.

Analysis of feature weights

After cross-validation and significance testing, we trained the model on all participants and extracted the feature weights for further analysis. First, we calculated the absolute value of the weights and summed them across all connections (edges) for a given region (node) to compare the overall contribution of each region to the model, irrespective of the sign of the feature weights (). Next, to evaluate the spatial pattern of the feature weights, we calculated the mean signed feature weight for each node, reflecting the directionality of the effect of the drug manipulation according to the trained model. We then used this feature map to assess the biological relevance of our trained model. Specifically, we calculated the spatial correlation between this pattern of nodal feature weights with two sets of cortical features. The first was the widely used principal gradient of macroscale cortical organization (), which places each cortical region on a continuum between unimodal (i.e., sensorimotor cortices) and transmodal (i.e., association cortices) function. The second set of cortical features was selected on the basis of the known pharmacology of benzodiazepines like alprazolam. Alprazolam is a positive allosteric modulator of the GABAA receptor, and of the six GABAA α subunits (α1 to α6), only subunits α1, α2, α3, and α5 are benzodiazepine sensitive (). To quantify the spatial distribution of the six GABAA α subunits, we extracted the microarray gene expression patterns for their corresponding GABAA receptor genes (GABRA1 to GABRA6) from the Allen Human Brain Atlas (data available at www.meduniwien.ac.at/neuroimaging/mRNA.html) (). For each of the six gene expression maps, we quantified the mean expression value within each cortical parcel and calculated the spatial correlation with the pattern of nodal SVM weights.

To test the significance of the spatial correlation between our pattern of cortical feature weights and each of the biological brain maps, we compared the observed correlation value to a null distribution generated with BrainSMASH (Brain Surrogate Maps with Autocorrelated Spatial Heterogeneity; https://brainsmash.readthedocs.io/) (). The spatial autocorrelation of brain maps can lead to inflated P values in spatial correlation analyses and must be accounted for in the creation of null models. BrainSMASH addresses this by generating permuted null brain maps that match the spatial autocorrelation properties of the input data. We used BrainSMASH to generate 10,000 spatial autocorrelation–preserving null permutations based on the input data and the pairwise distance matrix for the cortical parcellation, generating a null distribution of spatial correlation coefficients. We calculated two-tailed P values by squaring all correlation values (i.e., spatial R2) and calculating the proportion of times the null distribution exceeded the observed value.

Transmodal and unimodal classification models

Our primary hypothesis was that E:I ratio reductions would be specific to association cortices during youth. To test this hypothesis directly, we trained two additional models after applying an a priori feature selection step. Specifically, we thresholded the top and bottom quartiles of cortical parcels based on their position in the principal gradient of functional organization (), with the top 25% representing transmodal association cortex and the bottom 25% representing unimodal sensory cortex. We then created two new feature sets that restricted the input features to connections to these transmodal or unimodal areas only. This selection procedure ensured that the resulting numbers of features were equal between the two feature sets (9541 features per model). We then trained and validated the transmodal and unimodal models according to the procedures described above.

Developmental analyses

Application of the pharmacological model to the developmental dataset

After training and validating the pharmacological benzodiazepine models, we applied the models to the functional connectivity data for each participant in the developmental sample. For each participant, each model yielded the distance from the classification hyperplane that separates the two classes (drug versus placebo). Observations close to the hyperplane (distance values near zero) are less representative of the class, and those further from the hyperplane are more representative. The distance metric is such that values greater than zero indicate more drug-like patterns of functional connectivity, and values less than zero indicate more placebo-like patterns of connectivity. As the pharmacological effect of alprazolam is to increase GABAergic inhibitory signaling, more drug-like patterns reflect greater GABAergic inhibitory modulation of functional connectivity. As such, more drug-like patterns were interpreted to reflect a reduced E:I balance relative to more placebo-like patterns. These distance metrics were normally distributed and thus provided a continuous measure of E:I balance for use in further analyses. We first applied the model trained on all the input features and then applied the transmodal- and unimodal-specific models, generating three sets of distance values per participant.

Developmental regression models

To assess the developmental trajectory of E:I balance, we modeled the classification distance metrics from each model as a function of age using penalized splines within a GAM. GAMs allow us to flexibly capture linear or nonlinear age effects while penalizing overfitting. To test for windows of significant change across the age range, we calculated the first derivative of the smooth function of age from the GAM model using finite differences and then generated a simultaneous 95% confidence interval of the derivative using the gratia library () in R. Intervals of significant change were identified as areas where the simultaneous confidence interval of the derivative does not include zero. To test whether the effect of age on classification distance differed between the transmodal and unimodal SVM models, we calculated the residualized change () in transmodal versus unimodal distance scores by regressing the unimodal distance out of transmodal distance. We then regressed the residualized change score on age using a GAM. All models included sex as a covariate as well as head motion and attentiveness as covariates of no interest. Head motion was quantified as mean framewise displacement during the fMRI acquisition. Attentiveness was quantified as the number of response omissions during the emotion identification task; this covariate was included to control for potential effects of arousal on model performance as alprazolam can cause drowsiness. All GAMs were fit using the mgcv library () in R.

As previously described (, ), PNC participants underwent a clinical assessment of psychopathology. Multiple domains of psychopathology symptoms were evaluated using a structured screening interview (GOASSESS); we used these data to investigate whether dimensions of psychopathology moderated developmental reductions in E:I balance. As has been detailed in previous work (, ), factor scores were derived from the clinical assessments using a bifactor confirmatory factor analysis model that included a general factor for overall psychopathology as well as four specific factors that primarily represent anxious-misery (mood and anxiety) symptoms, psychosis spectrum symptoms, behavioral symptoms (conduct and attention-deficit hyperactivity disorder), and fear symptoms (phobias). All five factors are orthogonal and can be considered jointly in analysis of imaging data. To sample a broad range of psychopathology, we expanded our inclusion criteria to include individuals with a history of psychiatric hospitalization and those receiving pharmacological psychiatric treatment (N = 1018; ages 8 to 21.7; M = 15.0; SD = 3.23; male/female = 462/556). We analyzed these data in a GAM that included age-by-factor score interactions for each factor from the bifactor model. Interactions were fit as bivariate smooth interactions with penalized splines using tensor interaction smooths (“ti” in mgcv).