# Adolescent maturation of cortical excitation-inhibition ratio based on individualized biophysical network modeling

Amin Saberi, Kevin Wischnewski, Kyesam Jung, Leon Lotter, H. Lina Schaare, Tobias Banaschewski, Gareth Barker, Arun Bokde, Sylvane Desrivières, Herta Flor, et al.

► To cite this version:

Amin Saberi, Kevin Wischnewski, Kyesam Jung, Leon Lotter, H. Lina Schaare, et al. Adolescent maturation of cortical excitation-inhibition ratio based on individualized biophysical network modeling. Science Advances, 2025, 11 (23), pp.eadr8164. 10.1126/sciadv.adr8164. hal-05478387

HAL Id: hal-05478387

https://hal.science/hal-05478387v1

Submitted on 26 Jan 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

NEUROSCIENCE

# Adolescent maturation of cortical excitation-inhibition ratio based on individualized biophysical network modeling

Amin Saberi<fnref n="1" />,²,³, Kevin J. Wischnewski<fnref n="1" />,²,⁴, Kyesam Jung<fnref n="1" />,², Leon D. Lotter<fnref n="1" />,²,⁵, H. Lina Schaare<fnref n="1" />,²,³, Tobias Banaschewski⁶, Gareth J. Barker⁷, Arun L. W. Bokde⁸, Sylvane Desrivières⁹, Herta Flor¹⁰,¹¹, Antoine Grigis¹², Hugh Garavan¹³, Penny Gowland¹⁴, Andreas Heinz¹⁵,¹⁶, Rüdiger Brühl¹⁷, Jean-Luc Martinot¹⁸, Marie-Laure Paillère Martinot¹⁸,¹⁹, Eric Artiges¹⁸,²⁰, Frauke Nees⁶,¹⁰,²¹, Dimitri Papadopoulos Orfanos¹², Herve Lemaitre¹²,²², Luise Poustka²³, Sarah Hohmann⁶, Nathalie Holz⁶, Christian Baeuchl²⁴, Michael N. Smolka²⁴, Nilakshi Vaidya²⁵, Henrik Walter¹⁵,¹⁶, Robert Whelan²⁶, Gunter Schumann²⁵,²⁷,²⁸,²⁹, IMAGEN Consortium†, Tomáš Paus³⁰, Juergen Dukart<fnref n="1" />,², Boris C. Bernhardt³¹, Oleksandr V. Popovych<fnref n="1" />,², Simon B. Eickhoff<fnref n="1" />,², Sofie L. Valk<fnref n="1" />,²,³*

The excitation-inhibition ratio is a key functional property of cortical microcircuits which changes throughout an individual's lifespan. Adolescence is considered a critical period for maturation of excitation-inhibition ratio. This has primarily been observed in animal studies. However, there is limited human in vivo evidence for maturation of excitation-inhibition ratio at the individual level. Here, we developed an individualized in vivo marker of regional excitation-inhibition ratio in human adolescents, estimated using large-scale simulations of biophysical network models fitted to resting-state functional imaging data from both cross-sectional (n = 752) and longitudinal (n = 149) cohorts. In both datasets, we found a widespread decrease in excitation-inhibition ratio in association areas, paralleled by an increase or lack of change in sensorimotor areas. This developmental pattern was aligned with multiscale markers of sensorimotor-association differentiation. Although our main findings were robust across alternative modeling configurations, we observed local variations, highlighting the importance of methodological choices for future studies.

INTRODUCTION

The vast repertoire of cortical functions emerges from a careful tuning of the interactions between excitatory and inhibitory neurons in microcircuits embedded in the structural scaffolding of the brain (1). Excitation and inhibition, mainly transmitted via glutamate and γ-aminobutyric acid (GABA), respectively, are inseparable and

balanced, i.e., the inhibition generated in the cortical microcircuits is proportional to the local and incoming excitation (2). This phenomenon has been observed during both responses to external stimuli (3–5) and spontaneous cortical activity (4, 6). The excitation-inhibition (E-I) balance is proposed to be essential for central aspects of cortical functioning, including the dynamic stability of activity (7), efficient

Copyright © 2025 The
Authors, some rights
reserved; exclusive
licensee American
Association for the
Advancement of
Science. No claim to
original U.S.
Government Works.
Distributed under a
Creative Commons
Attribution License 4.0
(CC BY).

<footnote n="1">¹Institute of Neuroscience and Medicine - Brain and Behaviour (INM-7), Research Centre Jülich, Jülich, Germany. ²Institute of Systems Neuroscience, Medical Faculty and University Hospital Düsseldorf, Heinrich Heine University Düsseldorf, Düsseldorf, Germany. ³Otto Hahn Group Cognitive Neurogenetics, Max Planck Institute for Human Cognitive and Brain Sciences, Leipzig, Germany. ⁴Institute of Mathematics, Faculty of Mathematics and Natural Sciences, Heinrich Heine University Düsseldorf, Düsseldorf, Germany. ⁵Max Planck School of Cognition, Stephanstrasse 1A, 04103 Leipzig, Germany. ⁶Department of Child and Adolescent Psychiatry and Psychotherapy, Central Institute of Mental Health, Medical Faculty Mannheim, Heidelberg University, Square J5, 68159 Mannheim, Germany. ⁷Department of Neuroimaging, Institute of Psychiatry, Psychology & Neuroscience, King's College London, London, UK. ⁸Discipline of Psychiatry, School of Medicine and Trinity College Institute of Neuroscience, Trinity College Dublin, Dublin, Ireland. ⁹Social, Genetic and Developmental Psychiatry Centre, Institute of Psychiatry, Psychology & Neuroscience, King's College London, London, UK. ¹⁰Institute of Cognitive and Clinical Neuroscience, Central Institute of Mental Health, Medical Faculty Mannheim, Heidelberg University, Square J5, Mannheim, Germany. ¹¹Department of Psychology, School of Social Sciences, University of Mannheim, 68131 Mannheim, Germany. ¹²NeuroSpin, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France. ¹³Departments of Psychiatry and Psychology, University of Vermont, Burlington, VT 05405, USA. ¹⁴Sir Peter Mansfield Imaging Centre School of Physics and Astronomy, University of Nottingham, University Park, Nottingham, UK. ¹⁵Department of Psychiatry and Psychotherapy CCM, Charité – Universitätsmedizin Berlin, corporate member of Freie Universität Berlin, Humboldt-Universität zu Berlin and Berlin Institute of Health, Berlin, Germany. ¹⁶German Center for Mental Health (DZPG), site Berlin-Potsdam, Berlin, Germany. ¹⁷Physikalisch-Technische Bundesanstalt (PTB), Braunschweig and Berlin, Germany. ¹⁸Institut National de la Santé et de la Recherche Médicale, INSERM U1299 "Trajectoires développementales en psychiatrie", Université Paris-Saclay, Ecole Normale supérieure Paris-Saclay, CNRS, Centre Borelli, Gif-sur-Yvette, France. ¹⁹Department of Child and Adolescent Psychiatry, Pitié-Salpêtrière Hospital, AP-HP Sorbonne Université, Paris, France. ²⁰Psychiatry Department, EPS Barthélémy Durand, Etampes, France. ²¹Institute of Medical Psychology and Medical Sociology, University Medical Center Schleswig-Holstein, Kiel University, Kiel, Germany. ²²Institut des Maladies Neurodégénératives, UMR 5293, CNRS, CEA, Université de Bordeaux, 33076 Bordeaux, France. ²³Department of Child and Adolescent Psychiatry, Center for Psychosocial Medicine, University Hospital Heidelberg, Heidelberg, Germany. ²⁴Department of Psychiatry and Psychotherapy, Technische Universität Dresden, Dresden, Germany. ²⁵Centre for Population Neuroscience and Stratified Medicine (PONS), Department of Psychiatry and Psychotherapy, Charité Universitätsmedizin Berlin, Berlin, Germany. ²⁶School of Psychology and Global Brain Health Institute, Trinity College Dublin, Dublin, Ireland. ²⁷Centre for Population Neuroscience and Precision Medicine (PONS), Institute for Science and Technology of Brain-inspired Intelligence (ISTBI), Fudan University, Shanghai, China. ²⁸Department of Psychiatry, University of Cambridge, Cambridge, UK. ²⁹German Centre for Mental Health, Berlin, Germany. ³⁰Departments of Psychiatry and Neuroscience, Faculty of Medicine and Centre Hospitalier Universitaire Sainte-Justine, University of Montreal, Montreal, Quebec, Canada. ³¹Multimodal Imaging and Connectome Analysis Laboratory, McConnell Brain Imaging Centre, Montreal Neurological Institute and Hospital, McGill University, Montreal, Canada.</footnote>

<footnote>*Corresponding author. Email: s.valk@fz-juelich.de</footnote>

<footnote>†IMAGEN Consortium authors and affiliations are listed in the Supplementary Materials.</footnote>
coding of the information (8), sharp tuning of sensory stimuli (2), and generation of synchronous cortical oscillations in gamma and beta ranges (2, 6, 9, 10). Conversely, disturbed E-I balance can lead to cortical circuit dysfunctioning and is hypothesized as a key pathophysiological mechanism under various neuropsychiatric conditions such as schizophrenia, autism spectrum disorder, and epilepsy (11–15).

Adolescence is a critical developmental period with substantial changes in the brain including maturation of the E-I ratio (16–18). During this period, several important changes occur in the architecture and function of excitatory and inhibitory neurons and synapses, which together are suggested to lead to a recalibration of the E-I ratio (16). For instance, postmortem histology of adolescent brains has shown a pruning of excitatory synapses within the prefrontal cortex in rats (19, 20), nonhuman primates (21, 22), and humans (23–25). In addition, postmortem transcriptomic studies of the prefrontal cortex in animals and humans have indicated marked changes in the expression of genes involved in inhibitory neurons and GABAergic signaling, including parvalbumin (9, 26, 27) and GABA type A (GABAA) receptor subunits (28–31). These transcriptomic changes are accompanied by the maturation of inhibitory function with stronger and shorter inhibitory postsynaptic currents, as observed in the prefrontal cortex of nonhuman primates (9, 30), overall indicating a relative increase in inhibitory synaptic transmission in this area (16, 32).

Currently available evidence on the in vivo maturation of the E-I ratio in humans is limited, as the invasive methods used in animal studies are not feasible in humans. However, in vivo proxies of the E-I ratio have been proposed, relying on its putative macroscale functional consequences captured in functional imaging (33, 34) and electrophysiology (35–37) or through biochemical quantification of glutamatergic or GABAergic neurotransmitters using magnetic resonance spectroscopy (38–40). Such approaches are informative but lack a certain level of mechanistic insight and detail that is observed with, for example, direct measurement of the excitatory and inhibitory input currents as done in animal research. Furthermore, studies on the development of the E-I ratio are often focused on selected areas, primarily within the prefrontal cortex, and the knowledge on the regional patterns of E-I ratio maturation across the whole cerebral cortex is limited. Biophysical network modeling (BNM) of the brain is a promising computational technique that can bridge different scales of investigation at a whole-cortical level. It provides a tool to noninvasively derive mechanistic inferences about a hidden brain feature at the microscale, such as the E-I ratio, based on the observed (empirical) in vivo data at the macroscale and has provided valuable insights into brain (dys)function (41–47). In this approach, the dynamic spontaneous activity of brain areas is simulated using biologically realistic models that are informed by, for instance, the blood oxygen level-dependent (BOLD) signal measured during resting-state functional magnetic resonance imaging (rs-fMRI) (48–51).

In this study, we aimed to investigate the in vivo maturation of the regional E-I ratio in adolescents at an individual level. To achieve this, we applied the BNM approach on two independent cross-sectional and longitudinal neuroimaging datasets from the Philadelphia Neuro-developmental Cohort (PNC) and the IMAGEN study (52, 53). We performed large-scale simulations of individualized BNMs (44, 54, 55), in which models were informed by structural connectivity (SC) and functional imaging data of each subject. The subject-level precision of these models allowed for mapping the estimated E-I ratio specifically in each individual using simulations that best represented their empirical data and furthermore enabled studying within-subject

maturation longitudinally. This extended a previous study that used the BNM approach to study E-I ratio development in the PNC dataset at the level of age groups (56). We demonstrated replicable effects across the two datasets, indicating cross-sectional and longitudinal age-related increases in relative inhibition in the association areas and no significant changes or relative increases in excitation in the sensorimotor areas. This pattern of the E-I ratio maturation was aligned with the proposed sensorimotor-association axis of the cortical neurodevelopment (18). Subsequently, given that the simulation results may be affected by various modeling and analytical choices (57–59) or might be confounded by the variability of underlying structural connectome, as well as the noise within the simulations and parameter optimization, we extensively assessed and demonstrated the robustness of our simulation-based findings against these nuisances. Last, we contrasted our marker of the E-I ratio with alternative, previously used BNM-based markers, highlighting methodological and conceptual considerations regarding their usage and interpretation.

## RESULTS

### Overview

We included 752 adolescents from the cross-sectional PNC dataset [409 female; mean age: 15.3 ± 2.4 (10 to 19) years] (52) and 149 participants from the longitudinal IMAGEN study (72 female; mean age: 14.4 ± 0.4 years at the baseline and 18.9 ± 0.5 years at follow-up) (53). Subject/session diffusion-weighted imaging (DWI) and rs-fMRI data were used to generate individual matrices of (i) structural connectome based on the density of white matter streamlines, (ii) functional connectivity (FC) matrix as the correlation of the BOLD signals, and (iii) functional connectivity dynamics (FCD) matrix as a measure of how the FC dynamically evolves through sliding windows of time during the scan across 100 cortical areas (60). Hereafter, we refer to FC and FCD matrices derived from the imaging data as empirical FC and FCD to distinguish them from simulated FC and FCD.

Next, we performed individualized BNM simulations and parameter optimizations for each subject/session to estimate their regional measures of the E-I ratio based on their in vivo imaging data (Fig. 1). We applied the reduced Wong-Wang model (61), which models each node as coupled excitatory and inhibitory neuronal pools, where the excitatory neuronal pools of different nodes are interconnected through the individual-specific SC. The model was controlled by global and regional free parameters, which were fit to the empirical resting-state functional data of the target subject/session using the covariance matrix adaptation-evolution strategy (CMA-ES) optimization algorithm (62–64). This involved running a maximum of 33,600 simulations per subject/session using an efficient implementation of BNM simulations on graphical processing units (GPUs; https://cubnm.readthedocs.io). The model parameters included a global parameter *G*, which scales the strength of interregional coupling, in addition to regional parameters $w_i^{EE}$, $w_i^{EI}$, and $w_i^{IE}$, which characterize the connectivity weights between excitatory and inhibitory neuronal pools within each node. Motivated by recent developments of this model (65, 66), we let $w_i^{EE}$ and $w_i^{EI}$ to vary across nodes, i.e., they were computed, independent of each other, through weighted combinations of six fixed biological maps that represent microstructural, functional, transcriptomic, and neurochemical heterogeneity of the human cerebral cortex. These maps were obtained from independent healthy adult samples and included average
<figure><img src="image_3.png" /><figcaption>**Fig. 1. Overview.** Individualized BNM simulation-optimization (A to C) was performed to derive the subject/session-specific regional measures of the E-I ratio, defined as time-averaged in silico input current to the excitatory neurons, $\langle I_i^E \rangle$ (D). The model consists of coupled excitatory and inhibitory neuronal pools in each node, where the excitatory neuronal pools of brain nodes are interconnected through the structural connectome of the given subject/session [(A), left]. The model is controlled by a global parameter G, which adjusts interregional coupling, in addition to regional parameters $w_i^{EE}$, $w_i^{EI}$, and $w_i^{IE}$, which characterize the connection weights between excitatory and inhibitory neuronal pools within each node. In each simulation, G, $w_i^{EE}$, and $w_i^{EI}$ are set by the optimizer, while $w_i^{IE}$ is determined on the basis of the FIC algorithm (B). The covariance matrix adaptation-evolution strategy was used to optimize model parameters given empirical data of a subject/session (C). The optimization goal was to maximize the goodness of fit by tuning 15 free parameters, including G, as well as bias and coefficient terms used to determine $w_i^{EE}$ and $w_i^{EI}$ based on six fixed biological cortical maps (fig. S1). The goodness of fit of each simulation to the empirical functional data [(A), right] was assessed as the correlation of FC matrices subtracted by their absolute mean difference and the KS distance of FCD matrices derived from the simulated or empirical BOLD signal. After completion of two optimization runs, the optimal simulation with the best goodness of fit to the empirical functional data of the target subject/session was selected (C). Last, the in silico input current to the excitatory neuron of each node $I_i^E$ was averaged across simulation time, resulting in an E-I ratio map for each subject/session (D).</figcaption></figure>

T1-weighted/T2-weighted ratio (T1w/T2w), average cortical thickness, principal gradient of FC (FC G1), principal axis of gene expression (Gene PC1), and average N-methyl-D-aspartate (NMDA) and GABA type A/Bz (GABA<sub>A/BZ</sub>) receptor positron emission tomography maps (fig. S1) (67–75). Furthermore, in each simulation, $w_i^{IE}$ was determined on the basis of an analytical-numerical feedback inhibition control (FIC) algorithm, which aimed to maintain the firing rate of excitatory neurons within a biologically plausible range of 3 Hz (61, 66). Then, from the optimal simulations of each subject/session, we extracted the in silico input current to the excitatory neuron of each node $I_i^E$, averaged across simulation time, which resulted in individual-specific $\langle I_i^E \rangle$ maps. The $\langle I_i^E \rangle$ values reflect in silico estimates of the regional E-I ratio, defined as the relative level of excitation compared to the relative level of inhibition exerted onto excitatory neurons, given that $I_i^E$ results from

the combination of excitatory input currents to each node (from itself and from the excitatory neurons of the other nodes through the SC) balanced by local inhibitory currents. Therefore, an increase in $\langle I_i^E \rangle$ can be interpreted as a relative increase in excitation or decrease in inhibition, i.e., an increase in E-I ratio, within a model region.

### Cross-sectional age-related variation of the E-I ratio

Then, we studied cross-sectional age-related variation of the E-I ratio during adolescence in the PNC dataset. The individualized optimal simulations of the PNC dataset showed a goodness of fit of 0.259 ± 0.101 to the empirical data (fig. S2). On the basis of these simulations, we found widespread significant age-related decreases in E-I ratio in association areas within the frontal, parietal, and temporal lobes, in contrast to its age-related increases in visual and
sensorimotor areas as well as the left posterior insula, controlled for goodness of fit, sex, and in-scanner rs-fMRI motion and adjusted for multiple comparisons at a false discovery rate (FDR) of 5% (Fig. 2A). The effect sizes across these regions, partial correlation of age and E-I ratio controlled for the confounds, ranged between -0.255 and 0.156. We then assessed the within-sample stability of the age effects across 100 subsamples of the data, each including half of the total sample with 376 subjects. The unthresholded age effects on the E-I ratio across all pairs of subsamples showed a mean correlation coefficient of 0.863 ± 0.047, indicating the high within-sample stability of the observed age effects (Fig. 2C). Of note, assessing maturational differences of the E-I ratio between males and females, we found no FDR-corrected significant age-by-sex interactions.

## Longitudinal changes of the E-I ratio

To extend and assess the replicability of our findings in the cross-sectional PNC study, we next investigated the longitudinal maturation of the E-I ratio in the independent IMAGEN dataset, including 149 participants assessed at the ages of 14 and 19 years. The individualized optimal simulations in the IMAGEN had mean goodness-of-fit values of 0.266 ± 0.102 at the baseline and 0.231 ± 0.113 at the follow-up session (fig. S3). Within these simulations, we found a significant longitudinal age-related decrease in E-I ratio in widespread association areas within the frontal, parietal, and temporal lobes and a significant increase in visual areas, controlled for goodness of fit, sex, in-scanner rs-fMRI motion, and site and adjusted for multiple comparisons at an FDR of 5% (Fig. 3A). The effect sizes across these regions, calculated as the standardized mean difference

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Cross-sectional effect of age on the E-I ratio during adolescence. (A) Effect of age on the E-I ratio, showing its significant age-related decrease (blue) and increase (red) during adolescence in the PNC dataset, after removing outliers and controlling for the goodness of fit, sex, and in-scanner rs-fMRI motion, corrected for multiple comparisons using FDR. (B) Unthresholded map of the effect of age on the E-I ratio. (C) Distribution of correlation coefficients between E-I ratio age effect maps of all pairs of subsamples across 100 half-split subsamples of the dataset.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Longitudinal effect of age on the E-I ratio during adolescence. (A) Longitudinal effect of age on the E-I ratio, showing its significant decrease (blue) and increase (red) through adolescence, using a mixed-effects model with random intercepts for each subject, after removing outliers and controlling for goodness of fit, sex, in-scanner rs-fMRI motion, and site, corrected for multiple comparisons using FDR. (B) Unthresholded effect of age on the E-I ratio. (C) Distribution of correlation coefficients between E-I ratio age effect maps of all pairs of subsamples across 100 half-split subsamples of the dataset. (D) Conjunction of regions showing significant decreases in E-I ratio associated with age in the PNC and IMAGEN datasets. (E) Spatial coalignment [Pearson correlation (r) or cosine similarity (cos)] of longitudinal effects of age on the E-I ratio in IMAGEN with a cross-sectional effect of age on the E-I ratio in PNC.</figcaption></figure>
of the baseline to follow-up session in E-I ratio controlled for the confounds, ranged between -0.299 and 0.229. We next assessed the within-sample stability of age effects using 100 subsamples of the IMAGEN data, each including half of the total sample with 74 subjects. This resulted in a mean correlation of *r* = 0.750 ± 0.108 between the E-I ratio age effect maps across all pairs of subsamples (Fig. 3C). Furthermore, we found no FDR-corrected significant age-by-sex interactions, yielding no evidence for sex differences in the maturation of the E-I ratio. Of note, because the quality of the tractograms in the baseline session of the IMAGEN dataset was lower, in these simulations, we used the SC of the follow-up session in the models of both sessions. However, in a subset of subjects with adequate quality of tractograms in both baseline and follow-up sessions (*n* = 110; 52 female), using models with session-specific SCs resulted in largely similar effects of age on the E-I ratio (*r* = 0.779, *P*<sub>spin</sub> < 0.001; cos = 0.841, *P*<sub>spin</sub> < 0.001; fig. S4).

We next assessed the similarity of cross-sectional and longitudinal age-related variation of the E-I ratio observed in the two datasets. Conjunction of regions with significant age effects on the E-I ratio in the PNC and IMAGEN datasets revealed 33 regions in the association cortices showing a significant decrease in E-I ratio, whereas no region showed a significant replicable age-related increase in E-I ratio (Fig. 3D). The mean E-I ratio across these association regions, after regressing out the effects of confounds, showed a correlation coefficient of *r* = -0.232 with age in the PNC (T = -6.64, P < 0.001) and a standardized mean difference of -0.312 between the sessions in IMAGEN (T = -4.17, P < 0.001; fig. S5). Furthermore, the unthresholded map of longitudinal effects of age on the E-I ratio in IMAGEN (Fig. 3B) was significantly coaligned (*r* = 0.630, *P*<sub>spin</sub> < 0.001; cos = 0.782, *P*<sub>spin</sub> < 0.001; Fig. 3E) with the map of cross-sectional effects of age on the E-I ratio observed in PNC (Fig. 2B). Therefore, overall, across the two datasets, we observed replicable cross-sectional and longitudinal effects, indicating a developmental decrease in E-I ratio in the association areas in contrast to a lack of significant changes or an increase in E-I ratio in sensorimotor areas.

**The neurodevelopmental pattern of the E-I ratio coaligns with the sensorimotor-association axis of cortical organization**

Having observed differential effects of age on the E-I ratio across cortical areas, we next sought to investigate the embedding of this spatial neurodevelopmental pattern across different domains of cortical organization as well as developmental transcriptomics. We first studied the spatial coalignment of the maps of E-I ratio maturation with a previously proposed sensorimotor-association axis of cortical neurodevelopment and the multimodal cortical features it was composed of (fig. S6 and table S1) (18). The maps of E-I maturation observed in both datasets were significantly (*P*<sub>spin</sub> < 0.05) correlated with the sensorimotor-association axis map (PNC: *r* = -0.617; IMAGEN: *r* = -0.607) as well as several of its components, notably including FC G1 (PNC: *r* = -0.691; IMAGEN: *r* = -0.641) and T1w/T2w (PNC: *r* = 0.437; IMAGEN: *r* = 0.548; Fig. 4A and fig. S7A). Next, comparing the maps of E-I maturation across seven canonical resting-state networks (76), in both datasets, we observed more negative age effects in the default mode, limbic and frontoparietal networks compared to the somatomotor and visual networks (Fig. 4B and fig. S7B). These findings indicated the coalignment of the E-I ratio maturational pattern with the sensorimotor-association

axis of the cortex with a higher age-related relative increase in inhibition toward the association areas.

Last, we performed developmental transcriptomics enrichment analysis of the E-I ratio maturation maps. Using partial least squares regression with the gene expression maps obtained from the Allen Human Brain Atlas (71, 73), we identified the top 500 genes expressed higher toward the negative and positive ends of the E-I ratio maturation maps. Next, we investigated the developmental enrichment of the two sets of genes using specific expression analysis of the BrainSpan dataset (77), comparing them against null sets of genes expressed in alignment with spun surrogates of the E-I ratio maturation maps (1000 permutations). We found the genes expressed toward the negative ends of the E-I ratio maturation maps to be enriched in later stages of development, significantly during neonatal to adolescence (PNC) or early childhood (IMAGEN) stages, in contrast to the genes expressed toward the positive ends of the E-I ratio maturation maps that were enriched in earlier fetal stages of development, although not significantly (Fig. 4D and fig. S7D).

**Sensitivity analyses**

Thus far, we observed consistent age effects in the adolescent maturation of the E-I ratio in two independent cross-sectional and longitudinal datasets across a sensorimotor-association axis by using simulations of individualized BNMs. However, these simulation-based findings may be sensitive to various modeling and analytical choices (57–59) as well as confounding effects of the underlying structural connectome or noise. Therefore, we next assessed the sensitivity of the E-I ratio and its age-related changes to such nuisances, including the effects of the interindividual variability of SC, modeling configurations, and the randomness within the optimizer and the simulations. To reduce the computational costs, we limited these analyses to a random subsample of 200 subjects from the PNC dataset. Hence, the effect of age on the E-I ratio in the “main run” was recalculated in this subsample for the comparisons with the alternative runs (Fig. 5A).

**Interindividual variability of structural connectome**

Using subject-specific SCs in the main analyses enabled modeling of brain function within an individualized structural scaffold, which better represents each subject. However, this potentially introduces interindividual variability of SCs as a source of variability in E-I ratio, particularly given that *I*<sub>*i*</sub><sup>E</sup> is directly related to the SC in the model Eq. 1. As a result, the associations of age with model-derived features may be confounded by age-related variation of trivial features of SC, such as the node-wise strength. However, the effect of age on the E-I ratio was robust to controlling for the node-wise SC strength as an additional confound (*r* = 0.960, *P*<sub>spin</sub> < 0.001; cos = 0.977, *P*<sub>spin</sub> < 0.001; fig. S8A). Furthermore, when SC variability was eliminated by using an identical template SC in the BNM simulations of all the subjects, the effect of age on the E-I ratio was coaligned with the effects observed in the main run (*r* = 0.485, *P*<sub>spin</sub> < 0.001; cos = 0.440, *P*<sub>spin</sub> < 0.001; Fig. 5B and fig. S8B) despite a poor average test-retest reliability of E-I ratio compared to the main run {intraclass correlation coefficient (ICC): 0.319 ± 0.165; range: [-0.179, 0.632]; fig. S8B}.

**Parcellation**

Using an alternative Schaefer parcellation with higher granularity of 200 nodes (60), the effect of age on the E-I ratio was largely consistent with the main run based on 100 nodes (*r* = 0.596, *P*<sub>spin</sub> < 0.001; cos = 0.615, *P*<sub>spin</sub> < 0.001; Fig. 5C and fig. S9A), yet we observed the
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Embedding of the E-I developmental pattern in the PNC dataset along the sensorimotor-association axis. (A) Spatial correlation of the E-I ratio maturation map in the PNC dataset with the maps of the sensorimotor-association cortical axis based on Sydnor et al. (18) (fig. S6). Colored diamonds show statistically significant ($P_{spin} < 0.05$) positive (red) and negative (blue) spatial correlations. (B) Distribution of the E-I ratio maturation map across the canonical resting-state networks (F = 13.85, $P_{spin} < 0.001$). Post hoc tests (Bonferroni-corrected) showed significantly more positive age effects in the visual (VIS) and somatomotor (SMN) compared to the limbic (LIM), frontoparietal (FPN), and default mode networks (DMN), in addition to more positive age effects in the dorsal attention network (DAN) compared to DMN. (C) Bottom: Mean expression of the top 500 genes associated with the E-I ratio maturation map, split into sets of negatively associated ($n = 187$, blue) and positively associated ($n = 313$, red) genes. Top: Specific expression analysis of the two sets of genes across developmental stages in the cortex. The y axis shows the negative log of FDR-corrected P values. Asterisks denote significantly enriched developmental stages compared to null genes based on spin surrogate maps (1000 permutations) and after FDR adjustment. SA, sensorimotor association; Evo., evolutionary; CMR, cerebral metabolic rate; Glu., glucose; CBF, cerebral blood flow; NeuroSynth PC1, principal component of NeuroSynth meta-analytical maps; LTC G1, principal gradient of laminar thickness covariance.</figcaption></figure>

poor average test-retest reliability of E-I ratio compared to the main run {ICC: 0.356 ± 0.186; range: [-0.075, 0.836]; fig. S9B}.

## Heterogeneity of regional parameters

We found consistent effects of age on the E-I ratio when alternative subsets of the six biological maps were used to determine the heterogeneity of regional parameters $w_i^{EE}$ and $w_i^{EI}$ (Fig. 5D and fig. S10): (i) using only T1w/T2w and FC G1 maps ($r = 0.776$, $P_{spin} < 0.001$; $\cos = 0.863$, $P_{spin} < 0.001$), (ii) using only NMDA and GABA<sub>A/BZ</sub> maps ($r = 0.331$, $P_{spin} < 0.001$; $\cos = 0.630$, $P_{spin} < 0.001$), and (iii) using T1w/T2w, FC G1, NMDA, and GABA<sub>A/BZ</sub> maps ($r = 0.853$, $P_{spin} < 0.001$; $\cos = 0.916$, $P_{spin} < 0.001$). Notably, the average ICC of E-I ratio compared to the main run was respectively 0.698 ± 0.108 {range: [0.183, 0.881]}, 0.530 ± 0.133 {range: [0.203, 0.795]}, and 0.819 ± 0.068 {range: [0.483, 0.931]}, indicating its moderate to good average test-retest reliability when using alternative sets of maps.

Next, we asked to what extent the spatial pattern of E-I ratio maturation is influenced by the spatial pattern of the underlying heterogeneity maps. To investigate this, we performed an alternative run by using a set of six null maps, generated by random spinning of the

six original maps on the cortical surface. In this null run, we observed weak and significant correlation of the E-I maturation map with two of the six null maps but none of the original maps. Conversely, the E-I maturation map of the main run based on the original maps was correlated more strongly with five of six original maps but only one of the null maps. This indicates that the choice of heterogeneity maps can to some extent influence the findings, although the nature of this influence is not trivial (fig. S11A).

This raises the question of whether one could use map-free alternative models that are free of this influence. To address this, we tested two map-free models including (i) a “homogeneous” model (three free parameters), which assumes the homogeneity of regional parameters and led to a significant drop of the goodness of fit by an average of -0.113 ± 0.062 (fig. S11B), and (ii) a “node-based heterogeneous” model (201 free parameters), which allows regional parameters of the nodes to vary independently as separate free parameters. Using the same number of simulations as the main run to fit this model, we observed a significant drop of the goodness of fit by an average of -0.059 ± 0.069 (fig. S11B), likely due to a reduced
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Sensitivity analyses. The unthresholded effect of age on the E-I ratio observed in a random subsample of the PNC dataset (n = 200) using the default configurations (A) compared to age effects observed using alternative configurations, including the following: (B) using a fixed template SC based on the MICs dataset, (C) definition of nodes based on a Schaefer parcellation with higher granularity of 200 nodes, (D) using alternative subsets of biological maps to determine the heterogeneity of regional parameters, (E) including the interhemispheric connections in the goodness of fit, and (F) using an alternative Gaussian noise seed. In (B) to (F), the statistics indicate spatial coalignment [Pearson correlation (r) or cosine similarity (cos)] of each map with the E-I ratio maturation map observed using the default configurations (A). (G) Pooled partial correlation of age with the E-I ratio (controlled for goodness of fit, sex, and in-scanner rs-fMRI motion) across (A) to (F) based on random-effects meta-analyses.</figcaption></figure>

rate of optimizer convergence (fig. S11C). This indicates the need for a much larger number of simulations to fit this complex model, which would be infeasible in our current, subject-level BNM approach. Thus, we argue that although these map-free alternatives can avoid the influence of maps, they either are too simplistic and lack enough biological detail (“homogeneous”) or are too complex and difficult to fit at this scale (“node-based heterogeneous”), both leading to a worse model fit. This in turn highlights the necessity of an intermediate feasible solution which is to induce the constrained heterogeneity of the regional parameters based on a set of biological maps. The effect of age on the E-I ratio based on these map-free models is reported in fig. S11D for the interested reader, although we refrain from interpreting these findings given the aforementioned issues of these models.

**Inclusion of interhemispheric connections in the cost function**

When interhemispheric connections were considered in the cost function, the effect of age on the E-I ratio was consistent with the main run (r = 0.734, *P*<sub>spin</sub> < 0.001; cos = 0.532, *P*<sub>spin</sub> < 0.001; Fig. 5E and fig. S12), yet the E-I ratio showed poor average test-retest reliability compared to the main run {ICC: 0.484 ± 0.128; range: [0.071, 0.743]}.

**Conduction velocity**

Next, we assessed the potential impact of interregional conduction delay at the individual level by repeating the optimal simulations while adding conduction delay between model regions informed by the subject-specific tractograms. We calculated the ICC of E-I ratio between optimal simulations with and without conduction delay and found its good test-retest reliability with a mean of 0.977 ± 0.017 using a velocity of 1 m/s to 0.997 ± 0.002 using a velocity of 6 m/s (fig. S13).

**Optimization random seed**

Within the parameter space, multiple local optima could exist that feature different E-I ratio values. To assess this, we calculated the node-wise ICC of E-I ratio across the optimal points obtained from the two CMA-ES runs of each subject. We found a mean ICC of 0.946 ± 0.021 {range: [0.874 to 0.982]} across nodes, indicating the good test-retest reliability of E-I ratio across alternative optima (fig. S14). Furthermore, in a ground truth recovery analysis, we fitted the model to synthetic functional data with known parameters and regional E-I ratio values. The recovered E-I ratio values showed a Pearson correlation of 0.977 (*P* < 0.001) and an ICC of 0.948 with the ground truth (fig. S15). Together, these findings indicate a low
risk of optimizers' convergence to different local optima with respect to the E-I ratio.

### Gaussian noise random sequence

We next investigated the impact of the random sequence (seed) of Gaussian noise introduced into the simulations by performing three analyses:

1) Ground truth recovery with different noise seeds: When using a different noise random seed than the one used in the ground truth simulation, the recovery fits decreased (0.811 ± 0.018, across 50 alternative seeds), with a lower Pearson correlation (0.672 ± 0.132) and ICC (0.339 ± 0.150) between the recovered and ground truth E-I ratio values, compared to when the same noise seed was used (fig. S15).

2) Test-retest reliability of E-I ratio in the optimal simulations across different noise seeds: We repeated the optimal simulations of each subject obtained in the main run and based on the default noise seed, with 50 alternative noise seeds. The median node-wise ICC of E-I ratio between the simulations using the original versus alternative noise seeds had an average of 0.810 ± 0.100 {range: [0.368 to 0.929]} across regions (fig. S16).

3) Comparing optimization runs performed using the default versus alternative noise seed: Running the full simulation-optimization process with a different Gaussian noise random sequence than the main run, we observed the poor average test-retest reliability of regional E-I ratio estimates {ICC: 0.471 ± 0.172; range: [-0.149, 0.771]}, indicating a notable effect of noise sequence on individual estimates. Nevertheless, the age-related E-I ratio effects remained largely consistent ($r = 0.537, P_{spin} < 0.001$; $\cos = 0.755, P_{spin} < 0.001$; Fig. 5F and fig. S17). Furthermore, when age effects were assessed on the basis of the average E-I ratio values across the two runs using the default and alternative random seeds (fig. S17D), they remained highly consistent with the effects observed in either run separately ($r = 0.844, P_{spin} < 0.001$; $\cos = 0.912, P_{spin} < 0.001$ for the main run; $r = 0.854, P_{spin} < 0.001$; $\cos = 0.930, P_{spin} < 0.001$ for the alternative run).

Together, these findings suggest that while the Gaussian noise random sequence influences simulation outcomes, the observed effects of age on the E-I ratio remain largely robust to this influence.

### Most consistent effects of age on the E-I ratio across alternative modeling configurations

Given the variability of the effects of age on the E-I ratio across the alternative models described above (Fig. 5, A to F), we next aimed to reconcile these findings and identify the most consistent effects of age on the E-I ratio across the different models. We performed random-effects meta-analyses, at the level of each parcel, to pool the partial correlation of age with the E-I ratio observed across the main and alternative configurations. We found a significant pooled age-related decrease in E-I ratio in association areas in frontal, parietal, and temporal lobes and its pooled increase in sensorimotor areas (Fig. 5G). The heterogeneity of the observed effects differed across the cortex and was higher in the somatomotor and parietal regions compared to visual and frontal regions (fig. S18A).

### Alternative BNM-based measures of the E-I ratio

In vivo estimation of the E-I ratio based on BNMs has been the aim of several previous studies using this or similar models (47, 56, 78–80). Yet, there has been no consensus on the BNM-based measures of the E-I ratio and various measures have been proposed and used across studies. Here, we present our findings regarding alternative

BNM-based measures of the E-I ratio used in the previous literature and highlight some considerations regarding their usage.

First, the optimal model parameters have been commonly used as measures of the E-I ratio (47, 78–80). Variation of these parameters can be interpreted as a shift of the ratio toward higher excitation (e.g., increase in $G$ or $w^{EE}$) or higher inhibition (e.g., increase in $w^{EI}$ or $w^{IE}$). However, in a multidimensional model in which these parameters can simultaneously covary and may degenerate, the interpretation of their variations is not straightforward. Across optimal simulations of subjects in the PNC dataset, we found significant associations between optimal parameters, such as a negative association of $w^{EI}$ and $w^{IE}$, indicating that lower excitatory-to-inhibitory connection weights are accompanied by (compensatory) higher inhibitory-to-excitatory connection weights (Fig. 6A). These associations were also reflected in the effects of age on the parameters (Fig. 6B). For example, there was an inverse correlation between the unthresholded effects of age on $w_i^{EI}$ and $w_i^{IE}$ ($r = -0.706, P_{spin} < 0.001$). The observed covariance between model parameters and their age effects indicates that these age effects should not be interpreted in isolation, and on the basis of these data, the net effect of age on the E-I ratio remains ambiguous.

In our study, as a solution to this problem of degeneracy between model parameters, we focused on a state variable of model nodes within the optimal simulations, which, as a “final common pathway,” reflects the collective outcome of the various model parameters on the E-I ratio within each node. At a neuronal level, the E-I ratio is commonly defined as the ratio, or the balance, between excitatory and inhibitory currents, potentials, and conductance onto excitatory neurons (81, 82). Accordingly, in our study, we quantified the E-I ratio based on time-averaged input current onto the excitatory neurons, $\langle I_i^E \rangle$, which reflects the net difference of the excitatory and inhibitory currents onto these neurons. Of note, our measure differs from another BNM-based measure of the E-I ratio based on model state variables, which was used in a similar previous study (56): the ratio of the time-average excitatory synaptic gating variable, $\langle S_i^E \rangle$, to the time-average inhibitory synaptic gating variable, $\langle S_i^I \rangle$. By inspecting the time series of model state variables in an example optimal simulation, we found that $I_i^E(t)$ and $S_i^E(t)/S_i^I(t)$ are positively correlated ($R^2 = 0.469$). However, in subsequent analyses, we found notable differences between these two alternative BNM-based measures of the E-I ratio: (i) Assuming that the firing rate of excitatory neurons, $r^E$, is an outcome of the E-I ratio (61), which indicates low versus high states of activity (81), we expect a measure of the E-I ratio to positively correlate with it. Accordingly, in an example simulation and using an exponential generalized linear mixed-effects model, we found that $r^E(t)$ correlates positively across nodes and time with both $I_i^E(t)$ ($R^2 = 0.977$) and $S_i^E(t)/S_i^I(t)$ ($R^2 = 0.573$) but is more strongly correlated with $I_i^E(t)$ (fig. S19). This was expected given that model Eq. 3 directly relates $I_i^E(t)$ to $r_i^E(t)$. (ii) Next, given the optimal simulation of 40 randomly selected subjects of the PNC dataset, we performed perturbed simulations in which one of the model parameters was increased or decreased by 10%, pushing the simulation to an expected state of increased/decreased excitation/inhibition (e.g., the 10% increase in $G$ is expected to push the system toward higher excitation). We then used paired *t* tests to compare each measure of the E-I ratio before and after the perturbation (fig. S20) and found larger effects of perturbations on $\langle I_i^E \rangle$ (mean $|T| = 21.480 \pm 5.383$) compared to
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Optimal model parameter interrelation and association with age in the PNC dataset. (A) Left: Pearson correlation of model parameter G and brain-averaged values of regional parameters $w^{EE}$, $w^{EI}$, and $w^{IE}$ across subjects are shown. Asterisks denote statistically significant correlations. Right: Interrelation of regional values of parameters $w^{EE}$, $w^{EI}$, and $w^{IE}$ across nodes and subjects based on a linear mixed-effects model with random intercepts and slopes per each node. (B) Left: Effect of age on optimal parameter G. Points represent the residual of G for each subject after removing confounds. Right: Unthresholded effect of age on regional parameters $w^{EE}$, $w^{EI}$, and $w^{IE}$. Age in years (y)</figcaption></figure>

$\langle S_i^E \rangle / \langle S_i^I \rangle$ (mean $|T| = 6.994 \pm 2.726$). This shows that $\langle I_i^E \rangle$ may be more sensitive to capture the variations in the E-I ratio caused by parameter perturbations. (iii) Last, we studied the effects of age on $\langle S_i^E \rangle / \langle S_i^I \rangle$ in a random subsample of 200 subjects from the PNC dataset (same as data used in the “Sensitivity analyses” section) and, across different model configurations, found significant increases primarily in unimodal areas, while most of the association areas showed no changes or small decreases (fig. S21). Notably, the effect of age on $\langle S_i^E \rangle / \langle S_i^I \rangle$, compared to its effect on $\langle I_i^E \rangle$, was less sensitive to the choice of SC, parcellation, inclusion of interhemispheric connections, and random seeds in the simulation (mean ICC: 0.999) and optimization (mean ICC: 0.963) but more sensitive to the choice of maps.

These findings, combined with the commonly used definition of the E-I ratio in experimental and theoretical research at the neuronal level, suggest that $\langle I_i^E \rangle$ may be a more direct and interpretable measure of the E-I ratio using the BNM approach compared to the alternatives used in the literature. They also highlight how modeling choices and parameters can affect the outcomes of E-I ratio changes with age.

## DISCUSSION

In this study, we used large-scale simulations of biologically realistic and individualized BNMs to estimate regional E-I ratio based on in vivo imaging data and evaluated its maturation during adolescence. We found a developmental decrease in the E-I ratio (higher inhibition or lower excitation) in the association areas, while the sensorimotor areas showed a lack of significant changes or a developmental increase. This finding was supported by imaging data from two independent datasets and through investigating both cross-sectional, interindividual age-related variations of the E-I ratio, as well as its longitudinal, within-individual changes through

adolescence. Our observed pattern of regional variability in the E-I ratio development aligned with the sensorimotor-association axis of cortical organization and highlighted the divergence of early versus late developmental timing of the sensorimotor and association areas. We extensively tested the sensitivity of our findings to various modeling nuisances and choices and found that despite certain variations, the E-I ratio maturation pattern was largely robust to them. Last, we contrasted our simulation-based measure of the E-I ratio to the alternative measures used in the literature and highlighted important considerations on their interpretations.

We found a robust and replicable developmental pattern of decreased E-I ratio in the association areas, indicating a relative increase in inhibition or decrease in excitation. This observation is in line with several findings from previous animal and human studies. At the molecular level, transcriptomics and proteomics analyses have revealed periadolescence changes in the expression of genes related to excitation and inhibition, such as the NMDA receptor subunits (83); calcium-binding proteins parvalbumin, calretinin, and calbindin, which are expressed in different types of interneurons (9, 26, 27, 32); and GABAA receptor subunits (28–31). These molecular shifts mirror changes in neuronal functional properties. For instance, within the prefrontal cortex, there is an increase in the subunit composition of the GABAA receptors, from $\alpha_2$- to $\alpha_1$-containing receptors, which have a faster decay time, resulting in faster synaptic inhibition (31, 32, 84). Consistent with this, recording of pyramidal neurons of the prefrontal cortex in nonhuman primates has indicated an increase in the strength as well as shortening of the inhibitory postsynaptic currents (9, 30). At the same time, microscopic investigation of the pyramidal neurons in the prefrontal cortex across different species has revealed a marked reduction of the excitatory synaptic density during adolescence (19–23, 25, 85). In humans, a recent study on the PNC dataset estimated the E-I ratio in vivo by modeling multi-variate patterns of FC and assessing their (dis)similarity to the FC of
adults receiving alprazolam, a GABAergic agonist, and reported a significant developmental decrease in the E-I ratio, which was specific to association areas (33). Furthermore, magnetic resonance spectroscopy has been used in several human studies to assess the maturation of in vivo levels of glutamate and GABA, primarily in the frontal areas, yet has reported inconsistent findings (38–40, 86–89). Notably, this inconsistency has been attributed to the sensitivity of lower-field scanners to macromolecule contaminants, which undermines the reliability of the findings in earlier studies (39). However, a recent study used 7-T scanners and found an age-related decrease in glutamate despite stable levels of GABA in the dorsolateral prefrontal cortex and decreased levels of both glutamate and GABA, despite no changes in their ratio, in the anterior cingulate and anterior insular cortices (39).

In contrast to the widespread decrease in the E-I ratio we observed in association regions during adolescence, in sensorimotor areas, we found an increase or no significant age-related changes of the E-I ratio. Human cortical maturation is suggested to unfold across a sensorimotor-association axis, with a differential temporal patterning indicating earlier maturation of the sensorimotor areas in contrast to later and more protracted maturation of the association cortices (18). In line with this, we found that the spatial pattern of E-I ratio maturation across cortical areas coaligned with the proposed sensorimotor-association axis of neurodevelopment (18). We additionally indicated that the genes preferentially expressed in association areas (showing a maturational decrease in the E-I ratio) are more prominent in later stages of development. The sensorimotor-association neurodevelopmental variation has been observed in cortical maturation of macrostructural features (90), intracortical myelination (91, 92), white matter connectivity (93), and functional organization (18, 94, 95), parallel to the maturation of excitation and inhibition (18). Consistent with the findings in the prefrontal cortex (21, 22), accelerated pruning of excitatory synapses around puberty has been observed in sensorimotor areas as well (85, 96, 97), although synaptic pruning in association regions is protracted and peaks later than in sensorimotor areas (18, 24). In addition, the maturation of parvalbumin inhibitory interneurons in association areas is suggested to be more prolonged (18, 98). Given the differences in the neurodevelopmental timing along the sensorimotor-association axis, it may be that the E-I ratio matures earlier in the sensorimotor areas before adolescence, and hence, we did not find a maturational increase in relative inhibition in these areas during our study age period. In line with our observation, two other studies using human in vivo markers of E-I ratio reported a significant increase in inhibition markers in association areas despite no significant changes in sensorimotor regions (33, 40). Future research can further investigate the regional differences in the timing of E-I ratio maturation by extending our approach to a wider age range, including developmental stages before and after adolescence. The presumed hierarchical progression of E-I maturation from sensorimotor areas in early life to association areas in later stages of development is thought to have important functional consequences (99, 100). Specifically, the maturation of the inhibitory circuitry in an area is mirrored by a critical period of enhanced experience-dependent plasticity, which is thought to support a shift of activity from spontaneous to stimulus evoked and in turn enhance the signal-to-noise ratio in performing task-dependent computations of that region. The hierarchical progression of this maturational cascade from sensorimotor areas in early life to association areas in later stages of development is in turn suggested to support the maturation of lower-order

sensory and motor functions toward higher-order social and executive functions (99, 100). Myelination, which is thought to increase near the end of these critical periods of plasticity (99), follows a similar hierarchical maturational trajectory (18), and it will be intriguing for future studies to investigate the link between the maturation of E-I ratio and myelination.

There are important methodological considerations in using BNMs to estimate E-I ratio based on imaging data. Our study extends upon a recent modeling study, which found a widespread relative increase in inhibition across the cortex, most prominently in the sensorimotor areas, by using BNMs constructed for 29 age groups of the PNC dataset and based on a template SC of an adult sample (56). In contrast, here, we used large-scale simulations to construct individualized BNMs (44, 54, 55), which allowed more specific simulation-based mapping of the E-I ratio in each individual subject, and this also enabled studying changes of the E-I ratio longitudinally within the same individual. In addition, individualized BNMs are shown to enhance the reliability of model parameters and fingerprinting accuracy of the simulated data (55). While our findings within the association areas were in agreement with the study by Zhang *et al.* (56), indicating a maturational decrease in the E-I ratio, they diverged in the sensorimotor areas. We suspect that this divergence can be attributed to several differences of the two studies, which, in addition to the usage of group-level versus individualized BNMs, include different simulation-based markers of the E-I ratio as well as the methodological details of image processing, modeling, and optimization. However, we presented findings that highlighted that $I_i^E$, compared to $S_i^E/S_i^I$, as well as model parameters, may be a more direct BNM-based marker of the E-I ratio and be closer to its common definition in the literature as the ratio, or the balance, between excitatory and inhibitory currents, potentials, and conductance onto excitatory neurons (81, 82).

To further investigate the role of methodological choices and confounds in our findings, we performed extensive sensitivity analyses and showed that our key findings remain robust to such nuisances. By meta-analytically pooling the age effects on E-I ratio across different models, we found a developmental pattern consistent with our main findings. However, while the overall maturational patterns were stable, we observed notable variations in the simulation outcomes across these alternative models, highlighting the importance of modeling and analytical decisions when using BNMs to infer hidden features such as the E-I ratio, in particular at the level of individual data. One critical methodological consideration is the selection of heterogeneity maps (or their absence), which can influence the spatial pattern of the E-I ratio maturation, although in nontrivial ways. While map-free models might eliminate this potential confound, they may not be biologically valid or technologically feasible options, and we showed that they lead to worse fitting of the model to the empirical data. Specifically, a “homogeneous” map-free model assumes that the local microcircuits are uniform across the cortex, contradicting well-established regional differences, for example, in the distribution of excitatory and inhibitory neuronal subtypes (65, 101). Conversely, a “node-based heterogeneous” model is too complex and difficult to fit, requiring considerably more simulations per subject. This makes it infeasible to use such a complex model within our individualized modeling approach applied to many subjects and sessions. Given challenges of map-free models, map-based heterogeneity offers a pragmatic intermediate solution that enables constrained heterogeneity of the parameters while maintaining
its feasibility. Another key methodological consideration is the char-
acteristics of the Gaussian noise incorporated in BNMs to account
for the inherent stochasticity of the brain (102, 103). We found that
variations in the random sequence of noise influenced simulation
outcomes, yet the observed effects of age on the E-I ratio remained
largely robust. This highlights that future BNM studies should explic-
itly take the often-overlooked impacts of noise random sequence
into account, for example, by repeating the simulations across mul-
tiple random seeds to assess the stability of their findings.

Recognizing the influence of modeling choices on simulation outcomes, we emphasize the need for future research to systematically explore and characterize these effects. While dedicated methodological studies should continue investigating how different modeling decisions shape BNM-derived measures (57, 58), studies applying BNMs to specific neurobiological questions should also assess the impact of alternative modeling choices on their reported findings. However, comprehensive sensitivity analyses can be computationally demanding, making it impractical to explore all sources of variability across large datasets. To address this, future research should consider using more efficient simulation implementations, such as the GPU-based approach introduced in this study, or developing methods that estimate modeling influences without running numerical simulations, such as analytical approximations (66) or deep learning-based methods (104).

Given the methodological flexibility of BNM paradigms, establishing optimal methodological choices in mapping the E-I ratio using BNMs at the individual level remains an open challenge. Therefore, it will be crucial for the future research to investigate the following: (i) Comprehensively assessing how modeling choices may affect the simulation-derived E-I ratio measures. (ii) Developing approaches to enable more complex models while maintaining their feasibility and avoiding overfitting. Such models may include a “node-based heterogeneous” model or models that are more biologically detailed via incorporation of layer-wise neuronal subtypes and feedback/feedforward connections (105, 106), the modulatory neurotransmitter systems (46), or the dynamics of α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid receptors (107). (iii) Evaluating the replicability of our findings in alternative datasets, particularly with higher imaging quality and more extensive follow-ups, as well as testing the fitting of BNMs to additional empirical measures such as magneto/electroencephalography (44, 108). (iv) Experimentally validating the E-I ratio measures obtained through fitting of BNMs to the fMRI data against empirical measures obtained using electrical recordings or evaluating them in response to excitatory/inhibitory interventions.

Overall, this study provides in vivo evidence, based on individu-
alized BNMs, of a replicable and robust decrease in cortical E-I ratio
in association areas during adolescence. The normative maturation
of E-I ratio is suggested to have important functional consequences
(18, 99) and its dysmaturation is believed to be associated with
various neurodevelopmental disorders (11, 12). For example, the
neurodevelopmental model of schizophrenia suggests that aberrant
cortical maturation, particularly in the development of excitatory and
inhibitory functions, may contribute to the emergence of the disease
later in life (12, 13, 17, 109). Future studies should investigate the
clinical relevance of adolescent maturation of E-I ratio in relation
to the risk and diagnosis of neurodevelopmental disorders such as
schizophrenia, as this could potentially offer biomarkers for early
detection and intervention.

MATERIALS AND METHODS

This research complies with the ethical regulations as set by The In-
dependent Research Ethics Committee at the Medical Faculty of the
Heinrich Heine University Düsseldorf (study number 2018-317).
We used previously published data from sources that have received
ethics approval from their respective institutions (52, 53, 110–112).

Participants

We studied adolescents from two population-based datasets, in-
cluding the cross-sectional PNC (52, 110, 111) and the longitudinal
IMAGEN dataset (53). We selected subjects and follow-up sessions
within the age range of 10 to 19 years. In the IMAGEN cohort, with-
in this age range, phenotypic assessments were conducted at the
ages of 14, 16, and 19 and imaging data were acquired at the ages of
14 and 19. Next, we excluded subjects with poor quality of raw or
processed imaging data, as detailed below in the “Image processing
and quality control” section. In the IMAGEN, the subjects were ex-
cluded if the imaging data had poor quality in any of the baseline or
follow-up sessions, except for the DWI data for which subjects were
selected on the basis of the image quality of the follow-up session.
Our final sample consisted of 752 adolescents from the PNC (409
female; mean age: 15.3 ± 2.4 years) and 149 participants from IMAGEN
(72 female; baseline mean age: 14.4 ± 0.4 years; follow-up mean age:
18.9 ± 0.5 years). The PNC data were collected at a single center in
Philadelphia, while the IMAGEN data were acquired in five dif-
ferent centers across Europe, in Dresden (n = 58), Paris (n = 45),
Mannheim (n = 30), London (n = 14), and Dublin (n = 2). The PNC
included subjects with European American (n = 351), African
American (n = 332) and mixed/other (n = 69) ethnicities, and the
IMAGEN cohort consisted of subjects with Caucasian (n = 134),
mixed (n = 12), and non-Caucasian (n = 3) ethnicities.

Image acquisition

T1w, rs-fMRI, and DWI data were acquired using 3-T scanners from different manufacturers (PNC: Siemens; IMAGEN: Siemens, Philips, General Electric). For more details on the image acquisition parameters, we refer the reader to the respective publications for each dataset [PNC: table 1 in (52); IMAGEN: table S5 in (53)]. Of particular relevance to our study, the repetition and acquisition times of the rs-fMRI images were respectively 3 s and 6:18 min in the PNC and 2.2 s and 6:58 min in the IMAGEN.

**Image processing and quality control**

T1w structural magnetic resonance imaging

T1w MRI images were processed using the recon-all command of FreeSurfer (version 7.1.1; https://surfer.nmr.mgh.harvard.edu/), which includes brain extraction, registration to standard space, tissue segmentation, and cortical surface reconstruction (113, 114). The quality of the FreeSurfer output was controlled on the basis of the Euler characteristic, which represents the number of cortical surface defects before correction. We excluded outlier subjects with Euler characteristic greater than Q₃ + 1.5 × IQR (interquartile range) of their cohort (115). The FreeSurfer output and raw T1w images were subsequently used in the pipelines of rs-fMRI and DWI processing.

Resting-state functional magnetic resonance imaging

rs-fMRI images were preprocessed using fMRIprep (version 22.0.0; https://fmriprep.org/en/stable/), which performed brain extraction, image registration and motion correction, estimation of confounds,
and when possible, susceptibility distortion and slice time corrections (116). The latter two steps were omitted in the IMAGEN data given the unavailability of slice timing and variability of field-map formats across centers. We next further processed the output of fMRIprep, which involved the following steps: (i) applying Schaefer-100 parcellation (60) by taking the average signal of all vertices within each parcel at each time point, (ii) removing the first three volumes, (iii) high-pass temporal filtering >0.013 Hz, (iv) regressing out confounds including the average signal of white matter and cerebrospinal fluid voxels as well as 24 motion parameters (translation and rotation in the three directions, in addition to their squares, derivatives, and squares of derivatives) (117), and (v) scrubbing motion outliers, defined on the basis of root mean squared translation >0.25 mm. The scrubbing was done by setting the signal in motion outlier volumes to zero while Z-scoring the rest of the volumes. This approach, compared to discarding the motion outliers, preserves the temporal structure of the BOLD signal, which is important in calculating dynamic FC.

We excluded subjects/sessions with high rs-fMRI in-scanner motion defined as less than 4 min of scan remaining after scrubbing motion outlier volumes or a time-averaged root mean square >0.2 mm. In addition, we performed visual quality control of the fMRIprep output and excluded subjects with gross misregistration or incomplete field of view.

**Functional connectivity.** FC was calculated as the Pearson correlation of BOLD signal time series between the cortical areas.

**Functional connectivity dynamics.** The FCD matrix represents the temporal variability of the dynamic patterns of FC computed across sliding windows of time (118). To compute the FCD matrix, we initially calculated time-resolved FC matrices of sliding windows (PNC: size of 30 s and step of 6 s; IMAGEN: size of 30.8 s and step of 4.4 s). We discarded edge windows and windows with ≥50% motion outliers. Subsequently, we computed FCD as the correlation between lower triangular parts of window FC patterns. The distribution of values within the FCD matrix represents the amount of recurrence of time-resolved FC patterns.

### Diffusion-weighted imaging

DWI images were processed using Micapipe (version 0.1.1; https://micapipe.readthedocs.io/en/latest/) (119), which combines tools from FSL (version 6.0.0) (120) and MRTrix3 (version 3.0.0) (121). This pipeline performed DWI processing steps including rigid-body alignment of images, Marchenko-Pastur principal component analysis denoising, Gibbs ringing correction, motion and eddy current-induced distortion correction, nonuniformity bias correction, registration to the processed structural image, brain mask generation, and estimation of fiber orientation distributions using spherical deconvolution. Then, on each image, probabilistic tractography was performed using the iFOD2 algorithm to estimate 10 million streamlines (122). In addition, a track density image was computed using the iFOD1 algorithm with 1 million streamlines, which was used for the quality control (123). The quality control of the tractograms was done by visual inspection of the tractogram density images.

### Structural connectivity

The SC matrix for each subject was created using Micapipe by parcellating the tractogram using the Schaefer-100 parcellation map (60) non-linearly registered to the DWI space. We subsequently normalized each SC matrix by division by its mean × 100, resulting in an equal mean of 0.01 in all SCs.

In addition to subject- and session-specific SCs of the adolescents, we used higher-quality DWI data (3 T, three shells, 140 directions) of an adult sample of 50 healthy volunteers (MICs dataset, 23 female; mean age: 29.5 ± 5.6 years) to construct a template SC (112). To do so, we obtained the SC of individual MICs subjects preprocessed using Micapipe and calculated a group-averaged SC by taking the median of streamline counts in each edge. The template SC was subsequently normalized by its mean × 100, similar to the subject- and session-specific SCs of the adolescents.

In the individualized models reported in the main analyses, we used subject-specific SCs, but the adult template SC was used in the sensitivity analysis on the confounding effects of interindividual variability in SC. Of note, in IMAGEN, given the lower quality of tractograms in the baseline session, the main analyses were performed by using the follow-up SC of each subject for the modeling of functional data in both the baseline and follow-up sessions. However, in a supplementary analysis, we additionally performed simulations using session-specific SCs, within a subset of the IMAGEN subjects with adequate quality of tractograms in both sessions (n = 110, 52 female).

### Biophysical network modeling

Next, we performed BNM simulation-optimization at the level of each individual subject/session (Fig. 1).

### Model simulation

We simulated the spontaneous neuronal activity of 100 cortical regions from the Schaefer atlas (60) as network nodes regulated by the reduced Wong-Wang model and interconnected through the SC (61). In short, this model describes the activity of large ensembles of interconnected excitatory and inhibitory spiking neurons in each area by a dynamic mean field model as a reduced set of dynamic equations governing the activity of coupled excitatory and inhibitory pools of neurons. In this reduced model, the excitatory synaptic currents are mediated by the NMDA receptors and the inhibitory synaptic currents are mediated by the GABA<sub>A</sub> receptors. Within each cortical region, the excitatory and inhibitory neuronal pools are interconnected, and between regions, the excitatory neuronal pools are coupled through a scaled SC matrix.

**Model equations.** The model is mathematically described by a set of dynamic equations (61). The total input current (in nA) to each excitatory and inhibitory neuronal pool of each cortical node i, $I_i^{(E/I)}$, is calculated as

$$I_i^{E}(t) = W^{E}I_b + w_i^{EE}S_i^{E}(t) + GJ_{NMDA} \sum_{j=1}^{N} C_{ij}S_j^{E}(t) - w_i^{IE}S_i^{I}(t) \quad (1)$$

$$I_i^I(t) = W^I I_b + w_i^E S_i^E(t) - w^I S_i^I(t) \quad (2)$$

where $W^E I_b = 0.382$ nA and $W^I I_b = 0.267$ nA are the baseline input currents; $S_i^{(E/I)}$ denote the synaptic gating variables; n = 100 is the number of nodes; C is the SC matrix, which together with the NMDA receptor conductance, $J_{NMDA} = 0.15$ nA, and G (global coupling), a free parameter of the model, determines the excitatory input current transmitted from the other nodes; $w_i^{EE}$ is the recurrent excitatory connection weight; $w_i^{EI}$ indicates the excitatory-to-inhibitory connection weight; $w_i^{IE}$ is the inhibitory-to-excitatory connection weight; and $w^{II} = 1.0$ denotes recurrent inhibitory
connection weight. The local connection weights $w_i^{EE}$, $w_i^{EI}$, and $w_i^{IE}$ can vary across nodes and simulations through the free parameters of the model and the FIC, as described below.

The total input current received by each neuronal pool is subsequently transferred to $r_i^{(E/I)}$, firing rates in Hz, using the sigmoidal neural response function, $H^{(E/I)}$

$$r_i^E(t) = H^E(t) I_i^E = \frac{a^E I_i^E(t) - b^E}{1 - e^{-d^E [a^E I_i^E(t) - b^E]}} \quad (3)$$

$$r_i^I(t) = H^I(I_i^I) = \frac{a^I I_i^I(t) - b^I}{1 - e^{-d^I [a^I I_i^I(t) - b^I]}} \quad (4)$$

where $a^E = 310$ nC$^{-1}$ and $a^I = 615$ nC$^{-1}$ determine the slope of $H^{(E/I)}$; $b^E = 0.403$ nA and $b^I = 0.288$ nA define the thresholds above which the firing rates increase linearly with the input currents; and $d^E = 0.16$ and $d^I = 0.087$ determine the shape of $H^{(E/I)}$ curvature around $b^{(E/I)}$.

Last, the synaptic gating variables, $S_i^{(E/I)}$, follow

$$\frac{dS_i^E(t)}{dt} = -\frac{S_i^E(t)}{\tau_E} + [1 - S_i^E(t)]\gamma r_i^E(t) + \sigma\nu_i(t) \quad (5)$$

$$\frac{dS_i^I(t)}{dt} = -\frac{S_i^I(t)}{\tau_I} + r_i^I(t) + \sigma\nu_i(t) \quad (6)$$

where $S_i^E$ is mediated by NMDA receptors with a decay time constant $\tau_E = 0.1$ s and $\gamma = 0.641$, and $S_i^I$ is mediated by GABA receptors with a decay time constant $\tau_I = 0.01$ s. $\nu_i(t)$ is uncorrelated standard Gaussian noise with an amplitude of $\sigma = 0.01$ nA. $S_i^{(E/I)}$ is subsequently bound within the range [0, 1]. Note that the fixed parameters in Eqs. 1 to 6 are based on a previous paper by Deco et al. (61).

**Model free parameters.** The model is controlled by 15 free parameters, including G, as well as bias terms and coefficients for a fixed set of six biological maps, which together determine the regional values of $w_i^{EE}$ and $w_i^{EI}$. More specifically, $w_i^p, p \in \{EE, EI\}$ is calculated as

$$w_i^p = w_b^p \cdot \left( 1 + \sum_{k=1}^{6} c_k^p M_{ki} \right) \quad (7)$$

where $w_b^p$, the bias term, and $c_k^p$, a vector of six coefficients, are free parameters. M is a 6 by 100 matrix including the Z-scored biological maps in Schaefer-100 parcellation (fig. S1). These maps were based on independent samples of healthy individuals, obtained from neuromaps (74) and Hansen et al. (70) and included the following: (i) group-averaged T1w/T2w ratio map of the Human Connectome Project dataset (68, 69, 124), (ii) group-averaged cortical thickness map of the Human Connectome Project dataset, (iii) FC G1 (72), (iv) principal axis of Allen Human Brain Atlas gene expression data (Gene PC1) (71, 73), (v) NMDA receptor density positron emission tomography map (67, 70), and (vi) GABA$_{A/BZ}$ receptor density positron emission tomography map (70, 75). The resulting $w_i^{(EE/EI)}$ maps were subsequently shifted if needed to ensure $\min[w_i^{(EE/EI)}] \ge 0.001$.

We used the following ranges for the model free parameters: $G = [0.5, 4.0]$ and $w_b^{(EE/EI)} = [0.05, 0.75]$. The range of coefficients for each map was defined as $[-\frac{1}{\max(map)}, -\frac{1}{\min(map)}]$, corresponding to $c_{T1w/T2w} = [-0.48, 0.59]$, $c_{cortical thickness} = [-0.40, 0.39]$, $c_{FC G1} = [-0.59, 0.72]$, $c_{Gene PC1} = [-0.36, 0.48]$, $c_{NMDA} = [-0.49, 0.42]$, and $c_{GABAa/bz} = [-0.30, 0.32]$.

**Feedback inhibition control.** FIC was used to determine the regional values of $w_i^{IE}$ in each simulation, given the SC and other model parameters. The FIC algorithm aims to maintain a state of E-I balance in each region by adjusting $w_i^{IE}$ to satisfy an excitatory firing rate close to 3 Hz, which is suggested to be within the biological range (61). We used a two-stage implementation of the FIC by combining the original numerical implementation (61) with an analytical solution proposed by Demirtaş et al. (66). The latter solution analytically solves for $w_i^{IE}$ to satisfy the self-consistency of the model equations under the steady-state condition with $\langle r_i^E \rangle \approx 3$ Hz, corresponding to $\langle S^E \rangle \approx 0.164757$ and $\langle I^E \rangle \approx 0.37738$ nA

$$w_i^{IE} = \frac{W^E I_b + w_i^{EE} \langle S^E \rangle + G J_{NMDA} \langle S^E \rangle - \langle I^E \rangle}{\langle S_i^I \rangle} \quad (8)$$

in which the steady-state inhibitory synaptic gating variable $\langle S_i^I \rangle = H^I(\langle I_i^I \rangle)\tau_I$ was estimated by solving for $\langle I_i^I \rangle$ in

$$W^I I_b + w_i^{EI} \langle S^E \rangle - H^I(\langle I_i^I \rangle)\tau_I - \langle I_i^I \rangle = 0 \quad (9)$$

Subsequently, analytical estimates of $w_i^{IE}$ values were fed into the numerical implementation of FIC and were adjusted numerically (61). In this approach, given the analytical estimates of $w_i^{IE}$, the model Eqs. 1 to 6 are numerically integrated for a short period of 10 s and, subsequently, the average input current to the excitatory pool of each brain region, $\langle I_i^E \rangle$, is calculated. If $\langle I_i^E \rangle - \frac{b_E}{a_E}$ in a region exceeded the target -0.026 by more than 0.005 nA, $w_i^{IE}$ is up/down-regulated when the input current is higher/lower than the target, and the simulation is repeated with the adjusted $w_i^{IE}$ values in the next trial. This procedure is repeated for 10 trials or until the FIC target is satisfied in all nodes. Note that the maximum number of FIC numerical adjustment trials used here is lower than that of the original implementation to facilitate the scaling of the simulations. Furthermore, as the initial $w_i^{IE}$ values are estimated analytically rather than being fixed to 1 (as was done in the original implementation), a smaller number of trials will be needed.

**Hemodynamics model.** The simulated synaptic activity of the excitatory population of each node, $S_i^E$, was subsequently fed to the Balloon-Windkessel model of hemodynamics to simulate the BOLD signal (125). This model is mathematically described by the following system of differential equations with state variables $x$ (vasodilatory signal), $f$ (blood inflow), $v$ (blood volume), and $q$ (deoxyhemoglobin content)

$$\frac{dx(t)}{dt} = S^E(t) - \kappa x(t) - \gamma [f(t) - 1] \quad (10)$$

$$\frac{df(t)}{dt} = x(t) \quad (11)$$

$$\tau \frac{dv(t)}{dt} = f(t) - v_a^{\frac{1}{n}}(t) \quad (12)$$
$$ \tau \frac{dq(t)}{dt} = \frac{f(t)}{\rho} \left\{ 1 - \left[ (1-\rho)^{\frac{1}{f(t)}} \right] \right\} - \frac{q(t)v^{\frac{1}{\alpha}}(t)}{v(t)} \quad (13) $$

where $κ = \frac{1}{0.65} s^{-1}$ is the rate of signal decay, $γ = \frac{1}{0.41} s^{-1}$ is the rate of flow-dependent elimination, $τ = 0.98$ s is the hemodynamic transmit time, $α = 0.32$ is Grubb's exponent, and $ρ = 0.34$ is the resting oxygen extraction fraction. These parameters were based on a previous paper by Friston et al. (125). Last, on the basis of the model state variables, the BOLD signal is calculated as

$$ B(t) = V_0 \left\{ k_1[1-q(t)] + k_2 \left[ 1 - \frac{q(t)}{v(t)} \right] + k_3[1-v(t)] \right\} \quad (14) $$

in which $V_0 = 2\%$ is the resting blood volume fraction (125), and $k_1 = 3.72$, $k_2 = 0.527$, and $k_3 = 0.53$ are dimensionless parameters that were derived for 3-T scans (66, 126). Last, the simulated BOLD signal was downsampled to match the repetition time of the empirical rs-fMRI data, i.e., 3 s for the PNC and 2.2 s for IMAGEN.

**Numerical integration of the models.** For each simulation, the model equations were numerically integrated using the Euler method with a time step of 0.1 ms for the neuronal model (Eqs. 1 to 6) and a time step of 1 ms for the hemodynamic model (Eqs. 10 to 13). The model simulations were performed using in-house code (https://github.com/amnsbr/bnm promp and https://cubnm.readthedocs.io; see the “Data and materials availability” section for more details) on GPUs at JURECA-DC (127), Raven, or Juseless high performance/throughput computing systems. This GPU implementation enabled efficient parallelization of calculations for individual simulations (across GPU “blocks”) and the regions within each simulation (across GPU “threads”). To match the duration of empirical rs-fMRI scans, the simulations were done for a biological duration of 450 s, from which the first 30 s was discarded to ensure that the BNM system’s state has stabilized.

## Model evaluation

The goodness of fit of the simulated BOLD signal given a set of candidate parameters and SC matrix to a target empirical BOLD signal was evaluated on the basis of three measures of static and dynamic FC, following previous studies (56, 65):

**Static edge-level FC.** The simulated FC was calculated as the correlation of simulated BOLD signal time series between nodes. The correspondence of simulated and empirical FC patterns was evaluated by calculating the Pearson correlation coefficient between the lower triangles of the matrices ($FC_{corr}$), with larger values representing higher correspondence.

**Static global FC.** The absolute difference of the averaged simulated and empirical FC matrices across all the lower triangular edges ($FC_{diff}$) was calculated to assess the similarity of global FC strength, with smaller values showing higher correspondence.

**Dynamic FC.** The simulated FCD matrix was constructed by calculating the correlation of FC patterns between sliding windows of simulated BOLD signals, as described previously for the empirical data. The correspondence of simulated and empirical FCD distributions was calculated as the Kolmogorov-Smirnov (KS) distance of their lower triangular parts ($FCD_{KS}$), with smaller values showing higher similarity of the distributions.

Subsequently, these measures were combined into a single measure of goodness of fit as $FC_{corr} - FC_{diff} - FCD_{KS}$. Of note, in goodness-of-fit calculations, following Demirtaş et al. (66), we excluded the interhemispheric connections. However, we also performed a robustness analysis in which these connections were included in the goodness-of-fit calculations.

## Parameter optimization

The model's free parameters ($n = 15$) were fit to the empirical data of each subject/session using the CMA-ES optimization algorithm (62–64). CMA-ES is an efficient evolutionary optimization algorithm that features a set of $\Lambda$ particles exploring the parameter space collaboratively in an iterative process. The particles from each iteration, which are individual simulations with different free parameters, are regarded as a generation from which only the best particles are selected to form the descendant population of the next generation. Specifically, at each generation, the cost function of each particle is calculated, as described below. Then, a weighted mean of the best fitting [$\Lambda/2$] particles is calculated. Then, a new generation of particles is created by taking $\Lambda$ samples from a multivariate normal distribution centered around the weighted mean of the best fitting [$\Lambda/2$] particles from the previous generation. The covariance is determined by a matrix that is updated to take the location of the current best points into account. In this way, the search distribution is adapted iteratively toward a concentration around the optimal solutions. This iterative process is continued for a maximum of 80 generations, following Wischnewski et al. (64), and eventually, the optimal point across all generations is selected as the optimal parameters for the best fit of the simulation to the given target empirical data. We also applied an early termination rule in which the iteration was stopped if there was no improvement in the cost function greater than 0.005 over the past 30 generations.

The optimization goal was to maximize the goodness of fit while minimizing a penalty term. Particles were penalized if (i) the parameter of sampled particles fell outside the prespecified ranges, in which case the parameters were corrected and a penalty was added (63), or (ii) the $\langle r_i^E \rangle$ was outside the range of 2 to 4 Hz in any node, indicating the insufficiency of the FIC. For the latter, the FIC penalty was calculated as

$$ \text{FIC penalty} = \frac{2}{n} \sum_{i=1}^{n} \left( 1 - e^{-0.05 \left\langle r_i^E \right\rangle - 3} \right) \quad (15) $$

in which $n$ is the number of nodes, and the summation is done across nodes with out-of-range $\langle r_i^E \rangle$. Of note, we refrained from setting the success of FIC as a hard constraint, but through this penalty term, applied it as a soft constraint, to allow for interindividual variability of E-I balance while keeping it within a biologically viable range (61).

Given the relatively high dimensionality of the optimization problem in our model and to sufficiently cover the large parameter space, we chose $\Lambda = 210$ in the CMA-ES. Consequently, with maximum 80 generations, this involved performing a maximum of 16,800 simulations per run for each subject/session, which necessitated an efficient GPU-based implementation. In addition, for each given subject/session, the model simulation-optimization was repeated twice using different random seeds of the optimization to assess and reduce the risk of local optima. We then took the best fitting simulation across the two runs of each subject/session as their optimal simulation for the next step.

## Estimation of the E-I ratio in silico

Thus far, we described the procedure for deriving the optimal parameters that result in a simulation best fitting to the empirical rs-fMRI data of a given subject/session using individualized BNMs. Then, given the optimal simulation for each subject/session, we extracted an in silico measure of regional E-I ratio. To do so, we calculated the average of total input current to the excitatory neurons
of each region after discarding the initial 30 s of the simulation. This measure can be interpreted as an in silico marker of the regional E-I ratio, as $I_i^E$ (Eq. 1) results from the combination of excitatory input currents to each node (from itself and from the excitatory neurons of the other nodes through the SC) balanced by the inhibitory currents from the inhibitory neuron of the same node. Therefore, an increase in $\langle I_i^E \rangle$ can be interpreted as an increase in E-I ratio, i.e., a relative increase in excitation or decrease in inhibition. Furthermore, using a similar approach, we calculated $\langle S_i^E \rangle / \langle S_i^I \rangle$ as an alternative marker of the E-I ratio used in a previous study (56).

## Perturbed simulations

We performed a control analysis to assess the effect of known perturbations in the model parameters on the alternative markers of the E-I ratio, namely, $\langle I_i^E \rangle$ and $\langle S_i^E \rangle / \langle S_i^I \rangle$. In this analysis, we randomly selected 40 subjects of the PNC dataset and, for each subject, given their optimal simulations, we performed perturbed simulations in which one of the model parameters ($G$, $w^{EE}$, $w^{EI}$, or $w^{IE}$) was increased or decreased by 10%, while the other three parameters were fixed to the optimal values. Notably, in these simulations, when $w^{IE}$ was not perturbed, it was not readjusted using FIC and was fixed to the $w^{IE}$ values obtained by the FIC run on the original optimal simulation. Similarly, perturbation of $w^{IE}$ was done by a 10% increase or decrease of these original $w^{IE}$ values. This was to ensure that only one parameter in the model is perturbed and, therefore, the net direction of the effect of perturbation on the E-I ratio is predictable. Last, the effect of perturbation on the E-I ratio markers in each subject was assessed using paired *t* tests comparing the marker across nodes before and after the perturbation.

## Statistical analyses

### Age effects

Given the in silico regional measures of the E-I ratio for each subject/session, $\langle I_i^E \rangle$, we performed group-level univariate statistical analyses to investigate the effects of age on these measures. Linear regression models were used to study the effect of age on the E-I ratio, with the goodness of fit, sex, and rs-fMRI in-scanner motion (based on time-averaged root mean squared translation) as confounds. In IMAGEN, longitudinal variation of E-I measures across the two sessions was assessed using similar linear mixed-effects regression models with random intercepts per subject and inclusion of site as an additional confound. In each model, we excluded outliers with a dependent variable $\ge 3$ SDs above/below the mean. We adjusted for multiple comparisons across regions using FDR based on the Benjamini/Hochberg method ($q < 0.05$). Similar models were used to investigate the effects of age on optimal model parameters and $\langle S_i^E \rangle / \langle S_i^I \rangle$. We used statsmodels (https://www.statsmodels.org/stable/index.html) (128) to perform regressions and FDR adjustment.

### Within-sample stability of age effects using subsampling

In each dataset, we randomly selected 100 subsamples of the subjects (stratified by sex and age group; *n* within each subsample: 376 in PNC and 74 in IMAGEN) and investigated the effects of age on the E-I ratio separately in each subsample. Subsequently, we calculated the correlation of unthresholded age effect maps between all pairs of subsamples and reported its distribution as a measure of within-sample stability.

### Spatial association of maps

We evaluated the spatial association of an E-I ratio maturation map X with a target map Y (another E-I ratio maturation map or a brain

feature map) using Pearson correlation. In the comparison of two E-I ratio maturation maps, we additionally used cosine similarity given that Pearson correlation is insensitive to mean shifts and, thus, the changes in direction of the age effects. To account for spatial autocorrelation, the statistical significance of these associations was assessed by constructing nonparametric null distributions of Pearson correlation or cosine similarity calculated between X and spun surrogates of the target map Y. Spin permutation was implemented at the parcel level using the ENIGMA Toolbox (https://github.com/MICA-MNI/ENIGMA) (129) and based on 1000 permutations. This approach was used in assessing (i) the between-sample replicability of E-I ratio maturation between the PNC and IMAGEN datasets, (ii) the stability of E-I ratio maturation within a subsample of 200 subjects from the PNC dataset across the sensitivity analyses compared to the effects observed in the main run, and (iii) the spatial alignment of E-I ratio maturation maps (based on the main or alternative runs) with the maps determining the heterogeneity of regional parameters, as well as the sensorimotor-association axis of cortical organization described in a recent study by Sydnor et al. (18) and its components or their substitutes (some overlapping with the regional parameter heterogeneity maps) (fig. S6 and table S2) (68, 69, 71–73, 130–136).

### Distribution of the E-I ratio maturation maps across the canonical resting-state networks

We assessed the association of the E-I ratio maturation maps with the map of seven canonical resting-state networks (76) using one-way analysis of variance (ANOVA) with post hoc Bonferroni-corrected independent *t* tests. To control for spatial autocorrelation, we assessed the statistical significance of resulting *F* and *T* statistics using null distributions generated from 1000 spun surrogates based on parcel-level spinning implemented in the ENIGMA Toolbox (129).

### Partial least squares regression of gene expressions and their developmental enrichment

Regional microarray expression data were obtained from six post-mortem brains (one female; age: 24.0 to 57.0) provided by the Allen Human Brain Atlas (https://human.brain-map.org) (71). Data were processed with the abagen toolbox (https://abagen.readthedocs.io/en/stable/) (73) using the Schaefer-100 atlas (60). Gene expression data from the right hemisphere were excluded because of the sparsity of samples and a large number of regions with no expression data. We subsequently used scikit-learn (https://scikit-learn.org/stable/) (137) and performed partial least squares regression analysis to identify gene expression patterns with high spatial coalignment with the E-I ratio maturation maps within the left hemisphere. After selecting the top 500 genes with the highest absolute weights, we divided them into two sets of positively and negatively associated genes. Subsequently, using an online tool, we performed developmental specific expression analysis of these genes (http://dougherty-tools.wustl.edu/CSEAtool.html) (77). This tool uses Fisher's exact test to assess the overlap between the set of provided genes and pre-defined sets of genes, which are up-regulated in each developmental stage and brain structure, identified on the basis of the BrainSpan Atlas of the Developing Human Brain (http://www.brainspan.org). Here, for each set of genes, we reported the negative log of *P* values based on the specificity index (pSI) threshold of 0.05 within the cerebral cortex and FDR adjusted across the developmental stages. To assess the statistical significance of the observed developmental enrichment pattern, we constructed null distributions of negative
log values across developmental stages and positively and negatively associated genes by repeating the procedure described above for 1000 spun surrogate maps of E-I ratio maturation and performing developmental specific expression analysis on the resulting null sets of genes positively and negatively associated with the surrogate maps. The resulting spin-permutation *P* values were subsequently FDR adjusted.

### Test-retest reliability of E-I ratio across simulations

In the sensitivity analyses, we measured the node-level test-retest reliability of E-I ratio across simulations of the same subject between two simulation runs by measuring the median absolute deviation ICC of each region between the alternative simulations and across subjects. We interpreted ICC values as poor (<0.50), moderate (0.50 to 0.75), or good (≥0.75). To assess the effect of age on ICC measures, the ICC was additionally calculated separately in the younger and older subgroups of the PNC subsample subjects (split by median age), and the resulting ICC values across nodes were compared between the two age groups using paired *t* test.

### Association of optimal model parameters across subjects

The association of optimal regional parameters ($w^{EE}$, $w^{EI}$, and $w^{IE}$) across subjects and nodes was tested via linear mixed-effects regressions with random intercepts and slopes per each node. These regressions were performed using the lme4 R package (138). In addition, we used Pearson correlation to test the association of the optimal *G* with the mean of optimal regional parameters.

### Association of model state variables

We randomly selected one subject from the PNC dataset and evaluated the association of state variables within its optimal simulation across time points and nodes. The model state variables $r_i^E(t)$, $I_i^E(t)$, and $S_i^E(t)/S_i^I(t)$ were sampled every repetition time after the initial 30 s of the simulation was removed. Subsequently, as $S_i^E(t)/S_i^I(t)$ approaches infinity when $S_i^I(t)$ is close to zero, we excluded the data points at the top 2.5 percentile of $S_i^E(t)/S_i^I(t)$. Subsequently, we used mixed-effects models with or without a logarithmic linking function to test the linear or exponential associations of $I_i^E(t)$ with $r_i^E(t)$, $S_i^E(t)/S_i^I(t)$ with $r_i^E(t)$, and $S_i^E(t)/S_i^I(t)$ with $I_i^E(t)$ and reported the results of the model with a lower Akaike information criterion. The mixed-effects models included random intercepts and slopes and were implemented using the lme4 R package (138) and its Python interface in the pymer4 package.

### Random-effects meta-analysis of age effects across alternative modeling configurations

The sensitivity analyses described below resulted in several E-I ratio maturation maps based on the alternative configurations used for the BNM simulation-optimization runs. We performed random-effects meta-analyses, independently for each parcel and across seven main and alternative runs of map-based models (reported in Fig. 5, A to E), to calculate their pooled effects. The effect sizes used in these meta-analyses were partial correlations of age with the E-I ratio (after removing outliers and controlling for the goodness of fit, sex, and in-scanner rs-fMRI motion). The heterogeneity of effect sizes in each parcel was assessed using the *I*² index and Cochran Q test, where *P*(*Q*) < 0.05 indicates significant heterogeneity. In addition to the main meta-analyses across map-based models, we performed supplementary meta-analyses in which two map-free (“homogeneous” and “node-based heterogeneous”) models were additionally included. The PyMARE (https://pymare.readthedocs.io/en/latest/) package was used to perform all the meta-analyses.

### Sensitivity analyses

To assess the influences of modeling and analytical choices as well as the effects of interindividual variability of SC and noise, we performed a series of sensitivity analyses on a randomly selected subsample of 200 subjects (stratified by sex and age bin) from the PNC dataset and compared them to the main run.

### Interindividual variability of the structural connectome

We assessed the potential effect of the interindividual variability of SCs in our findings by performing the following analyses: (i) We studied the effect of age on the E-I ratio, additionally controlling for the SC strength of each node, calculated as the row-wise sum of the SC. (ii) We performed BNM simulation-optimization using subject-specific functional data as the target but with the template SC of the MICs dataset determining the connectivity of model nodes, thereby eliminating the interindividual variability of SCs as a potential source of variability in the regional E-I ratio. These higher-quality DWI data from an adult sample were chosen to additionally assess the robustness of our results to potential inaccuracies of subject-level SCs derived from relatively lower-quality DWI data of the adolescent datasets.

### Parcellation, heterogeneity of regional parameters, interhemispheric connections, and Gaussian noise seed

In these sensitivity analyses, we performed the BNM simulation-optimization using alternative modeling configurations, including (i) using the Schaefer-200 parcellation, (ii) using T1w/T2w and FC G1 as the heterogeneity maps, (iii) using NMDA and GABA<sub>A/BZ</sub> as the heterogeneity maps, (iv) using T1w/T2w, FC G1, NMDA, and GABA<sub>A/BZ</sub> as the heterogeneity maps, (v) using six null maps generated by randomly spinning the original maps together, as implemented in the ENIGMA Toolbox (126), (vi) assuming homogeneous regional parameters $w_i^{EE}$ and $w_i^{EI}$ in a three-parameter “homogeneous” model, (vii) assigning independent free parameters for each regional parameter of each node in a 201-parameter “node-independent” model, (viii) including interhemispheric connections in the goodness of fit and cost calculations, and (ix) using a different Gaussian noise seed than the default. Subsequently, we compared these optimal simulations derived from these alternative models with the main run in terms of goodness-of-fit measures, test-retest reliability of E-I ratio, and the effect of age on the E-I ratio. To assess the effect of Gaussian noise seed on the findings, we additionally calculated the effect of age on the E-I ratio based on average estimates from models using the default and alternative noise seeds. Furthermore, in the “node-based heterogeneous” model, we assessed the level of optimizer convergence based on the range (max – min) of goodness of fit across the particles in the last generation and compared it to that of the main run using paired *t* test.

### Effects of Gaussian noise seed and interregional conduction delay in optimal simulations of the main run

These analyses were performed by repeating the optimal simulation of each subject obtained in the main analyses with alternative configurations, including (i) using 50 different randomization seeds for generating the Gaussian noise injected into the model (Eqs. 5 and 6) and (ii) adding conduction delay in the signal transmission between the model nodes. For the latter, delay was calculated as the SC edge length obtained from the tractography of each subject, divided by a conduction velocity. For each subject, we performed six alternative simulations with variable conduction velocities in the range of {1, ..., 6} m/s. Of note, in simulations with conduction delay, a recent history of $S_i^E$ in all nodes needs to be stored in GPU
memory so that the input of node *j* to node *i* at time *t* can be determined on the basis of $S_j^E$ at delay$_{ij}$ time points ago. Accordingly, to reduce memory needed for storing this history, we performed these simulations by updating the global input to each node at intervals of 1 ms, instead of 0.1 ms used in the main analyses. Subsequently, we calculated the goodness-of-fit measures of the alternative simulations to the empirical data of each subject and compared them with the goodness-of-fit measures of the main simulation. In addition, we calculated the ICC of E-I ratio between the main simulation and each of the alternative simulations. In the case of Gaussian noise seeds, we combined the goodness of fit as well as ICC measures of the 50 different seeds by taking their median.

## Ground truth recovery analysis

We first generated synthetic functional data by running a ground truth simulation with known (and arbitrary) model parameters and using the template SC. Then, similar to the approach used for fitting the models to the real empirical data of subjects, we performed two BNM simulation-optimization runs aimed to maximize the model fit to the synthetic FC and FCD of the ground truth simulation. Next, we selected the best of the optima resulting from two runs as the “recovered optimal simulation” and compared it to the ground truth simulation in terms of goodness-of-fit measures and recovery of regional E-I ratio. For assessing the recovery of regional E-I ratio, we used Pearson correlation and ICC to compare the arrays of recovered versus ground truth simulations. We performed two alternative types of recovery optimization runs: (i) using the same Gaussian noise seed as the ground truth simulation or (ii) using 50 alternative Gaussian noise seeds that differed from the one used in the ground truth simulation.

## Supplementary Materials

This PDF file includes:

List of IMAGEN Consortium members

Figs. S1 to S21

Table S1

## REFERENCES AND NOTES

1. J. M. Shine, E. J. Müller, B. Munn, J. Cabral, R. J. Moran, M. Breakspear, Computational models link cellular mechanisms of neuromodulation to large-scale neural dynamics. *Nat. Neurosci.* **24**, 765–776 (2021).

2. J. S. Isaacson, M. Scanziani, How inhibition shapes cortical activity. *Neuron* **72**, 231–243 (2011).

3. J. S. Anderson, M. Carandini, D. Ferster, Orientation tuning of input conductance, excitation, and inhibition in cat primary visual cortex. *J. Neurophysiol.* **84**, 909–926 (2000).

4. M. Okun, I. Lampl, Instantaneous correlation of excitation and inhibition during ongoing and sensory-evoked activities. *Nat. Neurosci.* **11**, 535–537 (2008).

5. A. Y. Y. Tan, L. I. Zhang, M. M. Merzenich, C. E. Schreiner, Tone-evoked excitatory and inhibitory synaptic conductances of primary auditory cortex neurons. *J. Neurophysiol.* **92**, 630–643 (2004).

6. B. V. Atallah, M. Scanziani, Instantaneous modulation of gamma oscillation frequency by balancing excitation with inhibition. *Neuron* **62**, 566–577 (2009).

7. Y. K. Wu, C. Miehl, J. Gjorgjieva, Regulation of circuit organization and function through inhibitory synaptic plasticity. *Trends Neurosci.* **45**, 884–898 (2022).

8. S. Denève, C. K. Machens, Efficient codes and balanced networks. *Nat. Neurosci.* **19**, 375–382 (2016).

9. G. Gonzalez-Burgos, T. Miyamae, D. E. Pafundo, H. Yoshino, D. C. Rotaru, G. Hoftman, D. Datta, Y. Zhang, M. Hammond, A. R. Sampson, K. N. Fish, G. Bard Ermentrout, D. A. Lewis, Functional maturation of GABA synapses during postnatal development of the monkey dorsolateral prefrontal cortex. *Cereb. Cortex* **25**, 4076–4093 (2015).

10. V. S. Sohal, F. Zhang, O. Yizhar, K. Deisseroth, Parvalbumin neurons and gamma rhythms enhance cortical circuit performance. *Nature* **459**, 698–702 (2009).

11. J. H. Foss-Feig, B. D. Adkinson, J. L. Ji, G. Yang, V. H. Srihari, J. C. McPartland, J. H. Krystal, J. D. Murray, A. Anticevic, Searching for cross-diagnostic convergence: Neural mechanisms governing excitation and inhibition balance in schizophrenia and autism spectrum disorders. *Biol. Psychiatry* **81**, 848–861 (2011).

12. T. R. Insel, Rethinking schizophrenia. *Nature* **468**, 187–193 (2010).

13. E. T. Rolls, G. Deco, A computational neuroscience approach to schizophrenia and its onset. *Neurosci. Biobehav. Rev.* **35**, 1644–1653 (2011).

14. D. M. Treiman, GABAergic mechanisms in epilepsy. *Epilepsia* **42**, 8–12 (2001).

15. K. Xie, J. Royer, R. Rodriguez-Cruces, L. Horwood, A. Ngo, T. Arafat, H. Auer, E. Sahlas, J. Chen, Y. Zhou, S. L. Valk, S.-J. Hong, B. Frauscher, R. Pana, A. Bernasconi, N. Bernasconi, L. Concha, B. C. Bernhardt, Pharmaco-resistant temporal lobe epilepsy gradually perturbs the cortex-wide excitation-inhibition balance. *bioRxiv* 590555 [Preprint] (2024). https://doi.org/10.1101/2024.04.22.590555.

16. A. Caballero, A. Orozko, K. Y. Tseng, Developmental regulation of excitatory-inhibitory synaptic balance in the prefrontal cortex during adolescence. *Semin. Cell Dev. Biol.* **118**, 60–63 (2021).

17. G. D. Hoftman, D. Datta, D. A. Lewis, Layer 3 excitatory and inhibitory circuitry in the prefrontal cortex: Developmental trajectories and alterations in schizophrenia. *Biol. Psychiatry* **81**, 862–873 (2017).

18. V. J. Sydnor, B. Larsen, D. S. Bassett, A. Alexander-Bloch, D. A. Fair, C. Liston, A. P. Mackey, M. P. Milham, A. Pines, D. R. Roalf, J. Seidlitz, T. Xu, A. Raznahan, T. D. Satterthwaite, Neurodevelopment of the association cortices: Patterns, mechanisms, and implications for psychopathology. *Neuron* **109**, 2820–2846 (2021).

19. G. Gonzalez-Burgos, S. Kroener, A. V. Zaitsev, N. V. Povysheva, L. S. Krimer, G. Barrionuevo, D. A. Lewis, Functional maturation of excitatory synapses in layer 3 pyramidal neurons during postnatal development of the primate prefrontal cortex. *Cereb. Cortex* **18**, 626–637 (2008).

20. W. A. Koss, C. E. Belden, A. D. Hristov, J. M. Juraska, Dendritic remodeling in the adolescent medial prefrontal cortex and the basolateral amygdala of male and female rats. *Synapse* **68**, 61–72 (2014).

21. S. A. Anderson, J. D. Classay, F. Condé, J. S. Lund, D. A. Lewis, Synchronous development of pyramidal neuron dendritic spines and parvalbumin-immunoreactive chandelier neuron axon terminals in layer III of monkey prefrontal cortex. *Neuroscience* **67**, 7–22 (1995).

22. J. P. Bourgeois, P. S. Goldman-Rakic, P. Rakic, Synaptogenesis in the prefrontal cortex of rhesus monkeys. *Cereb. Cortex* **4**, 78–96 (1994).

23. P. R. Huttenlocher, Synaptic density in human frontal cortex — Developmental changes and effects of aging. *Brain Res.* **163**, 195–205 (1979).

24. P. R. Huttenlocher, A. S. Dabholkar, Regional differences in synaptogenesis in human cerebral cortex. *J. Comp. Neurol.* **387**, 167–178 (1997).

25. Z. Petanjek, M. Judaš, G. Šimić, M. R. Rašin, H. B. M. Uylings, P. Rakic, I. Kostović, Extraordinary neoteny of synaptic spines in the human prefrontal cortex. *Proc. Natl. Acad. Sci. U.S.A.* **108**, 13281–13286 (2011).

26. K. N. Fish, G. D. Hoftman, W. Sheikh, M. Kitchens, D. A. Lewis, Parvalbumin-containing chandelier and basket cell boutons have distinctive modes of maturation in monkey prefrontal cortex. *J. Neurosci.* **33**, 8352–8358 (2013).

27. S. J. Fung, M. J. Webster, S. Sivagnanasundaram, C. Duncan, M. Elashoff, C. S. Weickert, Expression of interneuron markers in the dorsolateral prefrontal cortex of the developing human and in schizophrenia. *Am. J. Psychiatry* **167**, 1479–1488 (2010).

28. D. Datta, D. Arion, D. A. Lewis, Developmental expression patterns of GABAA receptor subunits in layer 3 and 5 pyramidal cells of monkey prefrontal cortex. *Cereb. Cortex* **25**, 2295–2305 (2015).

29. S. G. Fillman, C. E. Duncan, M. J. Webster, M. Elashoff, C. S. Weickert, Developmental co-regulation of the beta and gamma GABAA receptor subunits with distinct alpha subunits in the human dorsolateral prefrontal cortex. *Int. J. Dev. Neurosci.* **28**, 513–519 (2010).

30. T. Hashimoto, Q. L. Nguyen, D. Rotaru, T. Keenan, D. Arion, M. Beneyto, G. Gonzalez-Burgos, D. A. Lewis, Protracted developmental trajectories of GABAA receptor α1 and α2 subunit expression in primate prefrontal cortex. *Biol. Psychiatry* **65**, 1015–1023 (2009).

31. G. D. Hoftman, D. A. Lewis, Postnatal developmental trajectories of neural circuits in the primate prefrontal cortex: Identifying sensitive periods for vulnerability to schizophrenia. *Schizophr. Bull.* **37**, 493–503 (2011).

32. A. Caballero, K. Y. Tseng, GABAergic function as a limiting factor for prefrontal maturation during adolescence. *Trends Neurosci.* **39**, 441–448 (2016).

33. B. Larsen, Z. Cui, A. Adebimpe, A. Pines, A. Alexander-Bloch, M. Bertolero, M. E. Calkins, R. E. Gur, R. C. Gur, A. S. Mahadevan, T. M. Moore, D. R. Roalf, J. Seidlitz, V. J. Sydnor, D. H. Wolf, T. D. Satterthwaite, A developmental reduction of the excitation-inhibition ratio in association cortex during adolescence. *Sci. Adv.* **8**, eabj8750 (2022).

34. S. Trakoshis, P. Martinez-Cañada, F. Rocchi, C. Canella, W. You, B. Chakrabarti, A. N. Ruigrok, E. T. Bullmore, J. Suckling, M. Markicevic, V. Zerbi, MRC AIMS Consortium, S. Baron-Cohen, A. Gozzi, M.-C. Lai, S. Panzeri, M. V. Lombardo, Intrinsic excitation-
inhibition imbalance affects medial prefrontal cortex differently in autistic men versus women. *elife* 9, e55684 (2020).

35. V. Medel, M. Irani, N. Crossley, T. Ossandón, G. Boncompte, Complexity and 1/f slope jointly reflect brain states. *Sci. Rep.* **13**, 21700 (2023).

36. T. M. Snijders, B. Milivojevic, C. Kemner, Atypical excitation-inhibition balance in autism captured by the gamma response to contextual modulation. *Neuroimage Clin.* **3**, 65–72 (2013).

37. P. J. Uhlhaas, F. Roux, E. Rodriguez, A. Rotarska-Jagiela, W. Singer, Neural synchrony and the development of cortical networks. *Trends Cogn. Sci.* **14**, 72–80 (2010).

38. C. Ghisleni, S. Bollmann, S.-S. Poil, D. Brandeis, E. Martin, L. Michels, R. L. O’Gorman, P. Klaver, Subcortical glutamate mediates the reduction of short-range functional connectivity with age in a developmental cohort. *J. Neurosci.* **35**, 8433–8441 (2015).

39. M. I. Perica, F. J. Calabro, B. Larsen, W. Foran, V. E. Yushmanov, H. Hetherington, B. Tervo-Clemmens, C.-H. Moon, B. Luna, Development of frontal GABA and glutamate supports excitation/inhibition balance from adolescence into adulthood. *Prog. Neurobiol.* **219**, 102370 (2022).

40. M. M. Silveri, J. T. Sneider, D. J. Crowley, M. J. Covell, D. Acharya, I. M. Rosso, J. E. Jensen, Frontal lobe γ-aminobutyric acid levels during adolescence: Associations with impulsivity and response inhibition. *Biol. Psychiatry* **74**, 296–304 (2013).

41. J. Cabral, H. M. Fernandes, T. J. Van Harteltvelt, A. C. James, M. L. Kringelbach, G. Deco, Structural connectivity in schizophrenia and its impact on the dynamics of spontaneous functional networks. *Chaos* **23**, 046111 (2013).

42. G. Deco, J. Cruzat, J. Cabral, G. M. Knudsen, R. L. Carhart-Harris, P. C. Whybrow, N. K. Logothetis, M. L. Kringelbach, Whole-brain multimodal neuroimaging model using serotonin receptor maps explains non-linear functional effects of LSD. *Curr. Biol.* **28**, 3065–3074.e6 (2018).

43. S. Froudist-Walsh, D. P. Bliss, X. Ding, L. Rapan, M. Niu, K. Knoblauch, K. Zilles, H. Kennedy, N. Palomero-Gallagher, X.-J. Wang, A dopamine gradient controls access to distributed working memory in the large-scale monkey cortex. *Neuron* **109**, 3500–3520.e13 (2021).

44. V. K. Jirsa, T. Proix, D. Perdikis, M. M. Woodman, H. Wang, J. Gonzalez-Martinez, C. Bernard, C. Bénar, M. Guye, P. Chauvel, F. Bartolomei, The Virtual Epileptic Patient: Individualized whole-brain models of epilepsy spread. *Neuroimage* **145**, 377–388 (2017).

45. K. Jung, E. Florin, K. R. Patil, J. Caspers, C. Rubbert, S. B. Eickhoff, O. V. Popovych, Whole-brain dynamical modelling for classification of Parkinson's disease. *Brain Commun.* **5**, fca331 (2022).

46. M. L. Kringelbach, J. Cruzat, J. Cabral, G. M. Knudsen, R. Carhart-Harris, P. C. Whybrow, N. K. Logothetis, G. Deco, Dynamic coupling of whole-brain neuronal and neurotransmitter systems. *Proc. Natl. Acad. Sci. U.S.A.* **117**, 9566–9576 (2020).

47. B. Park, S.-J. Hong, S. L. Valk, C. Paquola, O. Benkarim, A. R. I. Bethlehem, A. Di Martino, M. P. Milham, A. Gozzi, B. T. T. Yeo, J. Smallwood, B. C. Bernhardt, Differences in subcortico-cortical interactions identified from connectome and microcircuit models in autism. *Nat. Commun.* **12**, 2225 (2021).

48. M. Breakspear, Dynamic models of large-scale brain activity. *Nat. Neurosci.* **20**, 340–352 (2017).

49. G. Deco, M. L. Kringelbach, Great expectations: Using whole-brain computational connectomics for understanding neuropsychiatric disorders. *Neuron* **84**, 892–905 (2014).

50. O. V. Popovych, T. Manos, F. Hoffstaedter, S. B. Eickhoff, What can computational models contribute to neuroimaging data analytics? *Front. Syst. Neurosci.* **12**, 68 (2019).

51. P. Ritter, M. Schirner, A. R. McIntosh, V. K. Jirsa, The virtual brain integrates computational modeling and multimodal neuroimaging. *Brain Connect.* **3**, 121–145 (2013).

52. T. D. Satterthwaite, M. A. Elliott, K. Ruparel, J. Loughead, K. Prabhakaran, M. E. Calkins, R. Hopson, C. Jackson, K. Keefe, M. Riley, F. D. Mensh, P. Sleiman, R. Verma, C. Davatzikos, H. Hakonarson, R. C. Gur, R. E. Gur, Neuroimaging of the Philadelphia Neurodevelopmental Cohort. *Neuroimage* **86**, 544–553 (2014).

53. G. Schumann, E. Loth, T. Banaschewski, A. Barbot, G. Barker, C. Büchel, P. J. Conrod, J. W. Dalley, H. Flor, J. Gallinat, H. Garavan, A. Heinz, B. Ittlerman, M. Lathrop, C. Mallik, K. Mann, J.-L. Martinot, T. Paus, J.-B. Poline, T. W. Robbins, M. Rietschel, L. Reed, M. Smolka, R. Spanagel, C. Speiser, D. N. Stephens, A. Ströhle, M. Struve, the IMAGEN consortium, The IMAGEN study: Reinforcement-related behaviour in normal brain function and psychopathology. *Mol. Psychiatry* **15**, 1128–1139 (2010).

54. K. Bansal, J. Nakuci, S. F. Muldoon, Personalized brain network models for assessing structure–function relationships. *Curr. Opin. Neurobiol.* **52**, 42–47 (2018).

55. J. W. M. Domhof, S. B. Eickhoff, O. V. Popovych, Reliability and subject specificity of personalized whole-brain dynamical models. *Neuroimage* **257**, 119321 (2022).

56. S. Zhang, B. Larsen, V. J. Sydnor, T. Zeng, L. An, X. Yan, R. Kong, X. Kong, R. C. Gur, R. E. Gur, T. M. Moore, D. H. Wolf, A. J. Holmes, Y. Xie, J. H. Zhou, M. V. Fortier, A. P. Tan, P. Gluckman, Y. S. Chong, M. J. Meaney, G. Deco, T. D. Satterthwaite, B. T. T. Yeo, In vivo whole-cortex marker of excitation-inhibition ratio indexes cortical maturation and cognitive ability in youth. *Proc. Natl. Acad. Sci. U.S.A.* **121**, e2318641121 (2024).

57. J. W. M. Domhof, K. Jung, S. B. Eickhoff, O. V. Popovych, Parcellation-induced variation of empirical and simulated brain connectomes at group and subject levels. *Natw. Neurosci.* **5**, 798–830 (2021).

58. K. Jung, S. B. Eickhoff, O. V. Popovych, Tractography density affects whole-brain structural architecture and resting-state dynamical modeling. *Neuroimage* **237**, 118176 (2021).

59. O. V. Popovych, K. Jung, T. Manos, S. Diaz-Pier, F. Hoffstaedter, J. Schreiber, B. T. T. Yeo, S. B. Eickhoff, Inter-subject and inter-parcellation variability of resting-state whole-brain dynamical modeling. *Neuroimage* **236**, 118201 (2021).

60. A. Schaefer, R. Kong, E. M. Gordon, T. O. Laumann, X.-N. Zuo, A. J. Holmes, S. B. Eickhoff, B. T. T. Yeo, Local-global parcellation of the human cerebral cortex from intrinsic functional connectivity MRI. *Cereb. Cortex* **28**, 3095–3114 (2018).

61. G. Deco, A. Ponce-Alvarez, P. Hagmann, G. L. Romani, D. Mantini, M. Corbetta, How local excitation-inhibition ratio impacts the whole brain dynamics. *J. Neurosci.* **34**, 7886–7898 (2014).

62. N. Hansen, "Benchmarking a Bl-population CMA-ES on the BBOB-2009 function testbed," in *Proceedings of the 11th Annual Conference Companion on Genetic and Evolutionary Computation Conference: Late Breaking Papers* (ACM, 2009), pp. 2389–2396.

63. N. Hansen, The CMA Evolution Strategy: A tutorial. arXiv:1604.00772 (2016).

64. K. J. Wischnewski, S. B. Eickhoff, V. K. Jirsa, O. V. Popovych, Towards an efficient validation of dynamical whole-brain models. *Sci. Rep.* **12**, 4331 (2022).

65. G. Deco, M. L. Kringelbach, A. Arnatkeviciute, S. Oldham, K. Sabaroedin, N. C. Rogasch, K. M. Aquino, A. Fornito, Dynamical consequences of regional heterogeneity in the brain's transcriptional landscape. *Sci. Adv.* **7**, eabf4752 (2021).

66. M. Demirtas, J. B. Burt, M. Helmer, J. L. Ji, B. D. Adkinson, M. F. Glasser, D. C. Van Essen, S. N. Sotiriopoulos, A. Anticevic, J. D. Murray, Hierarchical heterogeneity across human cortex shapes large-scale neural dynamics. *Neurology* **101**, 1181–1194.e13 (2019).

67. M. Galovic, K. Erlandsson, T. D. Fryer, Y. T. Hong, R. Manavaki, H. Sari, S. Chetcuti, B. A. Thomas, M. Fisher, S. Septhon, R. Canales, J. J. Russell, K. Sander, E. Ärstad, F. I. Aigbirhoe, A. M. Groves, J. S. Duncan, K. Thielemans, B. F. Hutton, J. P. Coles, M. J. Koepp, NEST investigators, Validation of a combined image derived input function and venous sampling approach for the quantification of [$^{18}$F]GE-179 PET binding in the brain. *Neuroimage* **237**, 118194 (2021).

68. M. F. Glasser, S. N. Sotiriopoulos, J. A. Wilson, T. S. Coalson, B. Fischl, J. L. Andersson, J. Xu, S. Jbabdi, M. Webster, J. R. Polimeni, D. C. Van Essen, M. Jenkinson, WU-MHP Consortium, The minimal preprocessing pipelines for the Human Connectome Project. *Neuroimage* **80**, 105–124 (2013).

69. M. F. Glasser, D. C. Van Essen, Mapping human cortical areas in vivo based on myelin content as revealed by T1- and T2-weighted MRI. *J. Neurosci.* **31**, 11597–11616 (2011).

70. J. Y. Hansen, G. Shafiei, R. D. Markello, K. Smart, S. M. L. Cox, M. Nørgaard, V. Beliveau, Y. Wu, J.-D. Gallezot, E. Aumont, S. Servaes, S. G. Scala, J. M. DuBois, G. Wainstein, G. Bezgin, T. Funck, T. W. Schmitz, R. N. Spreng, M. Galovic, M. J. Koepp, J. S. Duncan, J. P. Coles, T. D. Fryer, F. I. Aigbirhoe, C. J. McGinnity, A. Hammers, J.-P. Soucy, S. Baillet, S. Guimond, J. Hietala, M.-A. Bedard, M. Leyton, E. Kobayashi, P. Rosa-Neto, M. Ganz, G. M. Knudsen, N. Palomero-Gallagher, J. M. Shine, R. E. Carson, L. Tuominen, N. Dagher, B. Misic, Mapping neurotransmitter systems to the structural and functional organization of the human neocortex. *Nat. Neurosci.* **25**, 1569–1581 (2022).

71. M. J. Hawrylycz, E. S. Lein, A. L. Guilozet-Bongaarts, E. H. Shen, L. Ng, J. A. Miller, L. N. van de Lagemaat, K. A. Smith, A. Ebbert, Z. L. Riley, C. Abajian, C. F. Beckmann, A. Bernard, D. Bertagnolli, A. F. Boe, P. M. Cartagena, M. M. Chakravarty, M. Chapin, J. Chong, R. A. Dalley, B. D. Daly, C. Dang, S. Datta, N. Dee, T. A. Dolbeare, V. Faber, D. Feng, D. R. Fowler, J. Goldy, B. W. Gregor, Z. Haradon, D. R. Haynor, J. G. Hohmann, S. Horvath, R. E. Howard, A. Jeromin, J. M. Jochim, M. Kinnunen, C. Lau, E. T. Lazarz, C. Lee, T. A. Lemon, L. Li, Y. Li, J. A. Morris, C. C. Overly, P. D. Parker, S. E. Parry, M. Reding, J. J. Royall, J. Schulkin, P. A. Sequeira, C. R. Slaughterbeck, S. C. Smith, A. J. Sodt, S. M. Sunkin, B. E. Swanson, M. P. Vawter, D. Williams, P. Wohnoutka, H. R. Zielke, D. H. Geschwind, P. R. Hof, S. M. Smith, K. Koch, S. G. N. Grant, A. R. Jones, An anatomically comprehensive atlas of the adult human brain transcriptome. *Nature* **489**, 391–399 (2012).

72. D. S. Margulies, S. S. Ghosh, A. Goulas, M. Falkiewicz, J. M. Huntenburg, G. Lang, G. Bezgin, S. B. Eickhoff, F. X. Castellanos, M. Petrides, E. Jefferies, J. Smallwood, Situating the default-mode network along a principal gradient of macroscale cortical organization. *Proc. Natl. Acad. Sci. U.S.A.* **113**, 12574–12579 (2016).

73. R. D. Markello, A. Arnatkeviciute, J.-B. Poline, B. D. Fulcher, A. Fornito, B. Misic, Standardizing workflows in imaging transcriptomics with the abagen toolbox. *elife* **10**, e72129 (2021).

74. R. D. Markello, J. Y. Hansen, Z.-Q. Liu, V. Bazinet, G. Shafiei, L. E. Suárez, N. Blostein, J. Seidlitz, S. Baillet, T. D. Satterthwaite, M. M. Chakravarty, A. Raznahan, B. Misic, neuromaps: Structural and functional interpretation of brain maps. *Nat. Methods* **19**, 1472–1479 (2022).

75. M. Nørgaard, V. Beliveau, M. Ganz, C. Svarer, L. H. Pinborg, S. H. Keller, P. S. Jensen, D. N. Greve, G. M. Knudsen, A high-resolution in vivo atlas of the human brain's benzodiazepine binding site of GABAA receptors. *Neuroimage* **232**, 117878 (2021).

76. B. T. T. Yeo, F. M. Krienen, J. Sepulcre, M. R. Sabuncu, D. Lashkari, M. Hollinshead, J. L. Roffman, J. W. Smoller, L. Zöllei, J. R. Polimeni, B. Fischl, H. L. R. Buckner, The organization of the human cerebral cortex estimated by intrinsic functional connectivity. *J. Neurophysiol.* **106**, 1125–1165 (2011).
77. J. D. Dougherty, E. F. Schmidt, M. Nakajima, N. Heintz, Analytical approaches to RNA profiling data for the identification of genes enriched in specific cells. *Nucleic Acids Res.* **38**, 4218–4230 (2010).

78. P. C. Klein, U. Ettinger, M. Schirner, P. Ritter, D. Rujescu, P. Falkai, N. Koutsouleris, L. Kambeitz-Ilkovic, J. Kambeitz, Brain network simulations indicate effects of neuregulin-1 genotype on excitation-inhibition balance in cortical dynamics. *Cereb. Cortex* **31**, 2013–2025 (2021).

79. A. Monteverdi, F. Palesi, A. Costa, P. Vitali, A. Pichieccio, M. Cotta Ramusino, S. Bernini, V. Jirsa, C. A. M. Gandini Wheeler-Kingshott, E. D'Angelo, Subject-specific features of excitation/inhibition profiles in neurodegenerative diseases. *Front. Aging Neurosci.* **14**, 868342 (2022).

80. G. J. Yang, J. D. Murray, X.-J. Wang, D. C. Glahn, G. D. Pearlson, G. Repovs, J. H. Krystal, A. Anticevic, Functional hierarchy underlies preferential connectivity disturbances in schizophrenia. *Proc. Natl. Acad. Sci. U.S.A.* **113**, E219–E228 (2016).

81. B. Haider, A. Duque, A. R. Hasenstaub, D. A. McCormick, Neocortical network activity in vivo is generated through a dynamic balance of excitation and inhibition. *J. Neurosci.* **26**, 4535–4545 (2006).

82. H. He, H. T. Cline, What is excitation/inhibition and how is it regulated? A case of the elephant and the wisemen. *J. Exp. Neurosci.* **13**, 1179069519859371 (2019).

83. M. A. Henson, A. C. Roberts, K. Salimi, S. Vadlamudi, R. M. Hamer, J. H. Gilmore, L. F. Jarskog, B. D. Philpot, Developmental regulation of the NMDA receptor subunits, NR3A and NR1, in human prefrontal cortex. *Cereb. Cortex* **18**, 2560–2573 (2008).

84. C. E. Duncan, M. J. Webster, D. A. Rothmond, S. Bahn, M. Elashoff, C. Shannon Weickert, Prefrontal GABAA receptor α-subunit expression in normal postnatal human development and schizophrenia. *J. Psychiatr. Res.* **44**, 673–681 (2010).

85. J. P. Bourgeois, P. Rakic, Changes of synaptic density in the primary visual cortex of the macaque monkey from fetal to adult stage. *J. Neurosci.* **13**, 2801–2820 (1993).

86. S. D. McKeon, M. I. Perica, A. C. Parr, F. J. Calabro, W. Foran, H. Hetherington, C.-H. Moon, B. Luna, Aperiodic EEG and 7T MRS evidence for maturation of E/I balance supporting the development of working memory through adolescence. *Dev. Cogn. Neurosci.* **6**, 101373 (2024).

87. S. Blüml, J. L. Wisnowski, M. D. Nelson Jr., L. Paquette, F. H. Gilles, H. C. Kinney, A. Panigrahy, Metabolic maturation of the human brain from birth through adolescence: Insights from in vivo magnetic resonance spectroscopy. *Cereb. Cortex* **23**, 2944–2955 (2013).

88. T. Gleich, R. C. Lorenz, L. Pöhland, D. Raufelder, L. Deserno, A. Beck, A. Heinz, S. Kühn, J. Gallinat, Frontal glutamate and reward processing in adolescence and adulthood. *Brain Struct. Funct.* **220**, 3087–3099 (2015).

89. M. Shimizu, Y. Suzuki, K. Yamada, S. Ueki, M. Watanabe, H. Igarashi, T. Nakada, Maturational decrease of glutamate in the human cerebral cortex from childhood to young adulthood: A 1H-MR spectroscopy study. *Diatr. Res.* **82**, 749–752 (2017).

90. R. A. I. Bethlehem, J. Seidlitz, S. R. White, J. W. Vogel, K. M. Anderson, C. Adamson, S. Adler, G. S. Alexopoulos, E. Anagnostou, A. Areces-Gonzalez, D. E. Astle, B. Auyeung, M. Ayub, J. Bae, G. Ball, S. Baron-Cohen, R. Beare, S. A. Bedford, V. Benegal, F. Beyer, J. Blangero, M. B. Cábez, J. P. Boardman, M. Borzage, J. F. Bayard, N. Bourke, V. D. Calhoun, M. M. Chakravarty, C. Chen, M. Chertavian, J. G. Chetelat, Y. S. Chong, J. H. Cole, A. Corvin, M. Costantino, E. Courchesne, F. Crivello, V. L. Cropley, J. Crosbie, N. Crossley, M. Delarue, R. Delorme, S. Desrivières, G. A. Devenyi, M. A. Di Biase, R. Dolan, K. A. Donald, G. Donohoe, K. Dunlop, A. D. Edwards, J. T. Elison, C. T. Ellis, J. A. Elman, L. Eyler, D. A. Fair, E. Feczko, P. C. Fletcher, P. Fonagy, C. E. Franz, G. Galan-Garcia, A. Gholipour, J. Giedd, J. H. Gilmore, D. C. Glahn, I. M. Goodyer, P. E. Grant, N. A. Groenewold, F. M. Gunning, R. E. Gur, R. C. Gur, C. F. Hammill, O. Hansson, T. Hedden, A. Heinz, R. N. Henson, K. Heuer, J. Hoare, B. Holla, A. J. Holmes, R. Holt, H. Huang, K. Im, J. Ipser, C. R. Jack, A. P. Jackowski, T. Jia, K. A. Johnson, P. B. Jones, D. T. Jones, R. S. Kahn, H. Karlsson, L. Karlsson, R. Kawashima, E. A. Kelley, S. Kern, K. W. Kim, M. G. Kitzbichler, W. S. Kremen, F. Lalonde, B. Landeau, S. Lee, J. Lerch, J. D. Lewis, J. Li, W. Liao, C. Liston, M. V. Lombardo, J. Lv, C. Lynch, T. T. Mallard, M. Marcelis, R. D. Markello, S. R. Mathias, B. Mazoyer, P. McGuire, M. J. Meaney, A. Mechelli, N. Medic, B. Misic, S. E. Morgan, D. Mothersill, J. Nigg, M. Q. W. Ong, C. Ortinau, R. Ossenkoppele, M. Ouyang, L. Palaniyappan, L. Paly, P. M. Pan, C. Pantelis, M. M. Park, T. Paus, Z. Pausova, D. Paz-Linares, A. P. Binette, K. Pierce, X. Qian, J. Qiu, A. Qiu, A. Raznahan, T. Rittman, A. Rodrigue, C. K. Rollins, R. Romero-Garcia, L. Ronan, M. D. Rosenberg, D. H. Rowitch, G. A. Salum, T. D. Satterthwaite, H. L. Schaare, R. J. Schachar, A. P. Schultz, G. Schumann, M. Schöll, D. Sharp, R. T. Shinohara, I. Skoog, C. D. Smyser, R. A. Sperling, D. J. Stein, A. Stolicyn, J. Suckling, G. Sullivan, Y. Taki, B. Thyreau, R. Toro, N. Traut, K. A. Tsvetanov, N. B. Turk-Browne, J. J. Tuulari, C. Tzourio, É. Vachon-Presreau, M. J. Valdes-Sosa, P. A. Valdes-Sosa, S. L. Valk, T. van Amelsvoort, S. N. Vandekar, L. Vasung, L. W. Victoria, S. Villeneuve, A. Villringer, P. E. Vértes, K. Wagstyl, Y. S. Wang, S. K. Warfield, V. Warrier, E. Westman, M. L. Westwater, H. C. Whalley, A. V. Witte, N. Yang, B. Yeo, H. Yun, A. Zalesky, H. J. Zar, D. Zettergren, J. H. Zhou, H. Ziauddeen, A. Zugman, X. N. Zuo, 3R-BRAIN, AIBL, Alzheimer's Disease Neuroimaging Initiative, Alzheimer's Disease Repository Without Borders Investigators, CALM Team, Cam-CAN,

CCNP, COBRE, cVEDA, ENIGMA Developmental Brain Age Working Group, Developing Human Connectome Project, FinnBrain, Harvard Aging Brain Study, IMAGEN, KNE96, The Mayo Clinic Study of Aging, NSPN, POND, The PREVENT-AD Research Group, VETSA, E. T. Bullmore, A. F. Alexander-Bloch, Brain charts for the human lifespan. *Nature* **604**, 525–533 (2022).

91. H. Grydeland, P. E. Vértes, F. Väša, R. Romero-Garcia, K. Whitaker, A. F. Alexander-Bloch, A. Bjørnerud, A. X. Patel, D. Sederevičius, C. K. Tarnes, L. T. Westlye, S. R. White, K. B. Walhov, A. M. Fjell, E. T. Bullmore, Waves of maturation and senescence in micro-structural MRI markers of human cortical myelination over the lifespan. *Cereb. Cortex* **29**, 1369–1381 (2019).

92. C. Paquola, R. A. Bethlehem, J. Seidlitz, K. Wagstyl, R. Romero-Garcia, K. J. Whitaker, R. Vos De Wael, G. B. Williams, P. E. Vértes, D. S. Margulies, B. Bernhardt, E. T. Bullmore, Shifts in myeloarchitecture characterise adolescent development of cortical gradients. *eLife* **8**, e50482 (2019).

93. D. J. Simmonds, M. N. Hallquist, M. Asato, B. Luna, Developmental stages and sex differences of white matter and behavioral development through adolescence: A longitudinal diffusion tensor imaging (DTI) study. *Neuroimage* **92**, 356–368 (2014).

94. D. S. Bassett, C. H. Xia, T. D. Satterthwaite, Understanding the emergence of neuropsychiatric disorders with network neuroscience. *Biol. Psychiatry Cogn. Neurosci. Neuroimaging* **3**, 742–753 (2018).

95. H.-M. Dong, D. S. Margulies, X.-N. Zuo, A. J. Holmes, Shifting gradients of macroscale cortical organization mark the transition from childhood to adolescence. *Proc. Natl. Acad. Sci. U.S.A.* **118**, e202448118 (2021).

96. N. Zecevic, J. P. Bourgeois, P. Rakic, Changes in synaptic density in motor cortex of rhesus monkey during fetal and postnatal life. *Brain Res. Dev. Brain Res.* **50**, 11–32 (1989).

97. N. Zecevic, P. Rakic, Synaptogenesis in monkey somatosensory cortex. *Cereb. Cortex* **1**, 510–523 (1991).

98. R. K. Reh, B. G. Dias, C. A. Nelson, D. Kaufer, J. F. Werker, B. Kolb, J. D. Levine, T. K. Hensch, Critical period regulation across multiple timescales. *Proc. Natl. Acad. Sci. U.S.A.* **117**, 23242–23251 (2020).

99. B. Larsen, B. Luna, Adolescence as a neurobiological critical period for the development of higher-order cognition. *Neurosci. Biobehav. Rev.* **94**, 179–195 (2018).

100. T. Toyoizumi, H. Miyamoto, Y. Yazaki-Sugiyama, N. Atapour, T. K. Hensch, K. D. Miller, A theory of the transition to critical period plasticity: Inhibition selectively suppresses spontaneous activity. *Neuron* **80**, 51–63 (2013).

101. X.-J. Wang, Macroscopic gradients of synaptic excitation and inhibition in the neocortex. *Nat. Rev. Neurosci.* **21**, 169–178 (2020).

102. G. Deco, V. Jirsa, A. R. McIntosh, O. Sporns, R. Kötter, Key role of coupling, delay, and noise in resting brain fluctuations. *Proc. Natl. Acad. Sci. U.S.A.* **106**, 10302–10307 (2009).

103. A. A. Faisal, L. P. J. Selen, D. M. Wolpert, Noise in the nervous system. *Nat. Rev. Neurosci.* **9**, 292–303 (2008).

104. J. D. Griffiths, Z. Wang, S. H. Ather, D. Momi, S. Rich, A. Diaconescu, A. R. McIntosh, K. Shen, Deep learning-based parameter estimation for neuropsychological models of neuroimaging data. bioRxiv 492664 [Preprint] (2022). https://doi.org/10.1101/2022.05.19.492664.

105. P. Wang, T. R. Knösche, A realistic neural mass model of the cortex with laminar-specific connections and synaptic plasticity – Evaluation with auditory habituation. *PLOS ONE* **8**, e77876 (2013).

106. A. F. Mejias, J. D. Murray, H. Kennedy, X.-J. Wang, Feedforward and feedback frequency-dependent interactions in a large-scale laminar network of the primate cortex. *Sci. Adv.* **2**, e1601335 (2016).

107. H. Sheheitli, V. Jirsa, Incorporating slow NMDA-type receptors with nonlinear voltage-dependent magnesium block in a next generation neural mass model: Derivation and dynamics. *J. Comput. Neurosci.* **52**, 207–222 (2024).

108. M. Schirner, A. R. McIntosh, V. Jirsa, G. Deco, P. Ritter, Inerring multi-scale neural mechanisms with brain network modelling. *eLife* **7**, e28927 (2018).

109. M. S. Keshavan, J. Giedd, J. Y. F. Lau, D. A. Lewis, T. Paus, Changes in the adolescent brain and the pathophysiology of psychotic disorders. *Lancet Psychiatry* **1**, 549–558 (2014).

110. M. E. Calkins, K. R. Merikangas, T. M. Moore, M. Burstein, M. A. Behr, T. D. Satterthwaite, K. Ruparel, D. H. Wolf, D. R. Roalf, F. D. Mentch, H. Qiu, R. Chiavacci, J. J. Connolly, P. M. A. Sleiman, R. C. Gur, H. Hakonarson, R. E. Gur, The Philadelphia Neurodevelopmental Cohort: Constructing a deep phenotyping collaborative. *J. Child Psychol. Psychiatry* **56**, 1356–1369 (2015).

111. M. E. Calkins, T. M. Moore, K. R. Merikangas, M. Burstein, T. D. Satterthwaite, W. B. Bilker, K. Ruparel, R. Chiavacci, D. H. Wolf, F. Mentch, H. Qiu, J. J. Connolly, P. A. Sleiman, H. Hakonarson, R. C. Gur, R. E. Gur, The psychosis spectrum in a young U.S. community sample: Findings from the Philadelphia Neurodevelopmental Cohort. *World Psychiatry* **13**, 296–305 (2014).

112. J. Royer, R. Rodríguez-Cruces, S. Tavakol, S. Larivière, P. Herholz, Q. Li, R. Vos de Wael, C. Paquola, O. Benkarim, B. Park, A. J. Lowe, D. Margulies, J. Smallwood, A. Bernasconi, N. Bernasconi, B. Frauscher, B. C. Bernhardt, An open MRI dataset for multiscale neuroscience. *Sci. Data* **9**, 569 (2022).
113. A. M. Dale, B. Fischl, M. I. Sereno, Cortical surface-based analysis. I. Segmentation and surface reconstruction. *Neuroimage* **9**, 179–194 (1999).

114. B. Fischl, M. I. Sereno, A. M. Dale, Cortical surface-based analysis. II: Inflation, flattening, and a surface-based coordinate system. *Neuroimage* **9**, 195–207 (1999).

115. J. Monereo-Sánchez, J. J. A. de Jong, G. S. Drenthen, M. Beran, W. H. Backes, C. D. A. Stehouwer, M. T. Schram, D. E. J. Linden, J. F. A. Jansen, Quality control strategies for brain MRI segmentation and parcellation: Practical approaches and recommendations - Insights from the Maastricht study. *Neuroimage* **237**, 118174 (2021).

116. O. Esteban, C. J. Markiewicz, R. W. Blair, C. A. Moodie, A. I. Isik, A. Erramuzpe, J. D. Kent, M. Goncalves, E. DuPre, M. Snyder, H. Oya, S. S. Ghosh, J. Wright, J. Durnez, R. A. Poldrack, K. J. Gorgolewski, fMRIPrep: A robust preprocessing pipeline for functional MRI. *Nat. Methods* **16**, 111–116 (2018).

117. K. J. Friston, S. Williams, R. Howard, R. S. Frackowiak, R. Turner, Movement-related effects in fMRI time-series. *Magn. Reson. Med.* **35**, 346–355 (1996).

118. E. C. A. Hansen, D. Battaglia, A. Spiegler, G. Deco, V. K. Jirsa, Functional connectivity dynamics: Modeling the switching behavior of the resting state. *Neuroimage* **105**, 525–535 (2015).

119. R. R. Cruces, J. Royer, P. Herholz, S. Larivière, R. Vos de Wael, C. Paquola, O. Benkarim, B. Park, J. Degré-Pelletier, M. C. Nelson, J. DeKraker, I. R. Leppert, C. Tardif, J.-B. Poline, L. Concha, B. C. Bernhardt, Micapipe: A pipeline for multimodal neuroimaging and connectome analysis. *Neuroimage* **263**, 119612 (2022).

120. M. Jenkinson, C. F. Beckmann, T. E. J. Behrens, M. W. Woolrich, S. M. Smith, FSL *Neuroimage* **62**, 782–790 (2012).

121. J.-D. Tournier, R. Smith, R. Raffelt, R. Tabbara, T. Dhollander, M. Pietsch, D. Christiaens, B. Jurrisen, C.-H. Yeh, A. Connelly, MRtrix:3: A fast, flexible and open software framework for medical image processing and visualisation. *Neuroimage* **202**, 116137 (2019).

122. J. D. Tournier, F. Calamante, A. Connelly, "Improved probabilistic streamlines tractography by 2nd order integration over fibre orientation distributions," in *Proceedings of the International Society for Magnetic Resonance in Medicine* (John Wiley & Sons Inc, 2010), vol. 1670.

123. J.-D. Tournier, F. Calamante, A. Connelly, MRtrix: Diffusion tractography in crossing fiber regions. *Int. J. Imaging Syst. Technol.* **22**, 53–66 (2012).

124. D. C. Van Essen, S. M. Smith, D. M. Barch, T. E. J. Behrens, E. Yacoub, K. Ugurbil, WU-Minn HCP Consortium, The WU-Minn Human Connectome Project: An overview. *Neuroimage* **80**, 62–79 (2013).

125. K. J. Friston, L. Harrison, W. Penny, Dynamic causal modelling. *Neuroimage* **19**, 1273–1302 (2003).

126. J. Heinzle, P. J. Koopmans, H. E. M. den Ouden, S. Raman, K. E. Stephan, A hemodynamic model for layered BOLD signals. *Neuroimage* **125**, 556–570 (2016).

127. P. Thörnig, JURECA: Data Centric and Booster Modules implementing the Modular Supercomputing Architecture at Jülich Supercomputing Centre. *J. Large-scale Res. Facilit.* **7**, A182 (2021).

128. S. Seabold, J. Perktold, "Statsmodels: Econometric and statistical modeling with python," in *Proceedings of the 9th Python in Science Conference* (SciPy, 2010), vol. 57, pp. 10–25080.

129. S. Larivière, C. Paquola, B. Park, J. Royer, Y. Wang, O. Benkarim, R. Vos de Wael, S. L. Valk, S. I. Thomopoulos, M. Kirschner, L. B. Lewis, A. C. Evans, S. M. Sisodiya, C. R. McDonald, P. M. Thompson, B. C. Bernhardt, The ENIGMA Toolbox: Multiscale neural contextualization of multisite neuroimaging data. *Nat. Methods* **18**, 698–700 (2021).

130. J. Hill, T. Inder, J. Neil, D. Dierker, J. Harwell, D. Van Essen, Similar patterns of cortical expansion during human development and evolution. *Proc. Natl. Acad. Sci. U.S.A.* **107**, 13135–13140 (2010).

131. P. K. Reardon, J. Seidlitz, S. Vandekar, S. Liu, R. Patel, M. T. M. Park, A. Alexander-Bloch, L. S. Clasen, J. D. Blumenthal, F. M. Lalonde, J. N. Giedd, R. C. Gur, R. E. Gur, J. P. Lerch, M. M. Chakravarty, T. D. Satterthwaite, R. T. Sinohara, A. Raznahan, Normative brain size variation and brain shape diversity in humans. *Science* **360**, 1222–1227 (2018).

132. S. N. Vaishnavi, A. G. Vlassenko, M. M. Rundle, A. Z. Snyder, M. A. Mintun, M. E. Raichle, Regional aerobic glycolysis in the human brain. *Proc. Natl. Acad. Sci. U.S.A.* **107**, 17757–17762 (2010).

133. T. D. Satterthwaite, R. T. Sinohara, D. H. Wolf, R. D. Hopson, M. A. Elliott, S. N. Vandekar, K. Rupale, M. E. Calkins, D. R. Roalf, E. D. Sennatgas, C. Jackson, G. Erus, K. Prabhakaran, C. Davatzikos, J. A. Detre, H. Hakonarson, R. C. Gur, R. E. Gur, Impact of puberty on the evolution of cerebral perfusion during adolescence. *Proc. Natl. Acad. Sci. U.S.A.* **111**, 8643–8648 (2014).

134. T. Yarkoni, R. A. Poldrack, T. E. Nichols, D. C. Van Essen, T. D. Wager, Large-scale automated synthesis of human functional neuroimaging data. *Nat. Methods* **8**, 665–670 (2011).

135. R. A. Poldrack, A. Kittur, D. Kalar, E. Miller, C. Seppa, Y. Gil, D. S. Parker, F. W. Sabb, R. M. Bilder, The cognitive atlas: Toward a knowledge foundation for cognitive neuroscience. *Front. Neuroinform.* **5**, 17 (2011).

136. A. Saberi, C. Paquola, K. Wagstyl, M. D. Hettwer, B. C. Bernhardt, S. B. Eichhoff, S. L. Valk, The regional variation of laminar thickness in the human isocortex is related to cortical hierarchy and interregional connectivity. *PLoS Biol.* **21**, e3002365 (2023).

137. F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, T. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, Scikit-learn: Machine learning in Python. *J. Mach. Learn. Res.* **12**, 2825–2830 (2011).

138. D. Bates, M. Mächler, B. Bolker, S. Walker, Fitting linear mixed-effects models using lme4. *J. Stat. Softw.* **67**, 1–48 (2015).

**Acknowledgments:** We gratefully acknowledge computing time on the supercomputer JURECA (124) at Forschungszentrum Jülich under grant no. 'eidev'. This research has been conducted using data from the dbGaP Study Accession: phs000607.v3.p2. Support for the collection of the data from the Philadelphia Neurogenomics Cohort (PNC) was provided by grant RC2MH089983 awarded to R. Gur and RC2MH089924 awarded to H. Hakonarson. Subjects were recruited and genotyped through the Center for Applied Genomics (CAG) at The Children's Hospital in Philadelphia (CHOP). Phenotypic data collection occurred at the CAG/CHOP and at the Brain Behavior Laboratory, University of Pennsylvania. **Funding:** A.S. and S.L.V. were funded by the Max Planck Society (Otto Hahn award) and Helmholtz Association's Initiative and Networking Fund under the Helmholtz International Lab grant agreement InterLabs-0015 and the Canada First Research Excellence Fund (CFREF Competition 2, 2015-2016) awarded to the Healthy Brains, Healthy Lives initiative at McGill University, through the Helmholtz International BigBrain Analytics and Learning Laboratory (HIBALL). L.D.L. was supported by the Federal Ministry of Education and Research (BMBF) and the Max Planck Society (MPG), Germany. S.B.E. was supported by the Deutsche Forschungsgemeinschaft (DFG, El 816/21-1), the National Institute of Mental Health (R01-MH074457), and the European Union's Horizon 2020 Research and Innovation Programme under grant agreement no. 945539 (HBP SGA3). The IMAGEN Consortium received support from the following sources: the European Union-funded FP6 Integrated Project IMAGEN (Reinforcement-related behaviour in normal brain function and psychopathology) (LSHM-CT-2007-037286 to G.S.); the Horizon 2020 funded ERC Advanced Grant 'STRATIFY' (Brain network based stratification of reinforcement-related disorders) (695313 to G.S.); Horizon Europe 'environmentMENTAL' (101057429); UK Research and Innovation (UKRI) Horizon Europe funding guarantee (10041392 and 10038599 to G.S.); Human Brain Project [HBP SGA 2 (785907) and HBP SGA 3 (945539) to G.S.]; The Chinese government via the Ministry of Science and Technology (MOST); The German Center for Mental Health (DZPG); The Bundesministerium für Bildung und Forschung (BMBF grants 01GS08152 and 01EV0711, Forschungsnetz AERIAL 01EE1406A and 01EE1406B, and Forschungsnetz IMAC-Mind 01GL1745B to G.S.); The Deutsche Forschungsgemeinschaft (DFG grants SM 80/7-2, SFB 940, TRR 265, and NE 1383/14-1 to G.S.); The Medical Research Foundation and Medical Research Council (grants MR/R00465X/1 and MR/S020306/1 to S.D.); the National Institutes of Health (NIH) funded ENIGMA grants 5US4E0B20403-05, 1R56AG058854-01, and 0US4E0B20403 to S.D., as well as NIH R01DA049238 to G.S.; the National Institutes of Health, Science Foundation Ireland (16/ERCD/3797 to R.W.); and NSFC (grant 82150710554). Further support to the IMAGEN Consortium was provided by grants from the following: the ANR (ANR-12-SAMA-0004 and AAPG2019-GeBra to J.-L.M.), the Etanet Neuron (AF12-NEUR0008-01-WM2NA and ANR-18-NEUR00002-01-ADORE to J.-L.M.), the Fondation de France (00081242 to J.-L.M.), the Fondation pour la Recherche Médicale (DPA20140629802 to J.-L.M.), the Mission Interministérielle de Lutte-contre-les-Drogues-et-les-Conduites-Addictives (MILDECA to J.-L.M.), the Assistance-Publique-Hôpitaux-de-Paris and INSERM (interface grant to M.-L.P.M.), Paris Sud University IDEX 2012 to J.-L.M., the Fondation de l'Avene (grant AP-RM-17-013 to M.-L.P.M.), and the Fédération pour la Recherche sur le Cerveau. **Author contributions:** Conceptualization: A.S. and S.L.V. Methodology: A.S., K.J.W., K.J., O.V.P., and S.L.V. Software: A.S. and K.J.W. Formal analysis: A.S. Investigation: A.S., D.L.D., T.B., G.J.B., A.L.W.B., S.D., H.F., AG, H.G., P.G., A.H., R.B., J.-L.M., M.-L.P.M., E.A., F.N., D.P.O., H.L., L.P., S.H., N.H., C.B., M.N.S., N.V., H.W., R.W., G.S., and T.P., and IMAGEN Consortium. Resources: S.B.E., S.L.V., and A.S. Data curation: A.S., D.L.D., T.B., G.J.B., A.L.W.B., S.D., H.F., AG, H.G., P.G., A.H., R.B., J.-L.M., M.-L.P.M., E.A., F.N., D.P.O., H.L., L.P., S.H., N.H., C.B., M.N.S., N.V., H.W., R.W., G.S., T.P., and IMAGEN Consortium. Writing—original draft: A.S. and S.L.V. Writing—review and editing: all authors. Visualization: A.S. Supervision: S.L.V. and S.B.E. Project administration: S.L.V. and S.B.E. Funding acquisition: S.B.E. and S.L.V.T.B., G.J.B., A.L.W.B., S.D., H.F., A.G., H.G., P.G., A.H., R.B., J.-L.M., M.-L.P.M., E.A., F.N., D.P.O., H.L., L.P., S.H., N.H., C.B., M.N.S., N.V., H.W., R.W., G.S., and T.P. were the principal investigators of the IMAGEN Consortium. The data in the IMAGEN project, including the methodology, funding, and supervision needed for collecting the data were provided by the IMAGEN Consortium. **Competing interests:** T.B. served in an advisory or consultancy role for eye level, Infectopharm, Medice, Neurim Pharmaceuticals, Oberberg GmbH, and Takeda. He received conference support or speaker's fee from Janssen, Medice, and Takeda. He received royalties from Hogrefe, Kohlhammer, CIP Medien, and Oxford University Press. The present work is unrelated to the above grants and relationships. G.J.B. has received honoraria from General Electric Healthcare for teaching on scanner programming courses. L.P. served in an advisory or consultancy role for Roche and Viforpharm and received speaker's fees from Shire. She received royalties from Hogrefe, Kohlhammer, and Schattauer. The present work is unrelated to the above grants and relationships. The other authors declare that they have no competing interests. **Data and materials availability:** All the code used in this study for image processing, BNM, and the statistical analyses is available in a GitHub repository (https://github.com/amnsfr/eidev), which is deposited at https://doi.org/10.5281/zenodo.12097801. Notably, while most of the BNMs in this study were performed using https://github.com/amnsfr/bnm_cuda (deposited at https://doi.org/10.5281/zenodo.12097805), we recommend and have included scripts to run the simulation-optimization using the cuBNM toolbox (https://github.com/amnsfr/cuBNM and https://cubnm.readthedocs.io; deposited at https://doi.org/10.5281/zenodo.14246606), which we developed using a similar core code but with an improved interface and performance. In
addition, the cuBNM toolbox was used to perform some of the simulations reported in the manuscript, including the Schaefer-200 model, "node-based heterogeneous" model, tests of delay, and simulation seed effects, as well as parameter perturbation and ground-truth recovery analyses. The PNC and IMAGEN data are publicly available in (i) the Database of Genotypes and Phenotypes (dbGaP), accession number phs000607.v3.p2 ("Neurodevelopmental Genomics: Trajectories of Complex Phenotypes"), and (ii) the IMAGEN database (https://imagen2.cea.fr). Access to these datasets is restricted to authorized users and can be obtained by application to (i) dbGaP (https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs000607.v3.p2) and (ii) the IMAGEN dataset (https://imagen-project.org/the-imagen-dataset/). The

individual-level raw and processed data of MICs dataset is openly available at https://osf.io/j532r/. All the summary statistics are available in the project repository at https://doi.org/10.5281/zenodo.12097081. All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials.

Submitted 17 July 2024

Accepted 25 April 2025

Published 4 June 2025

10.1126/sciadv.adr8164