Abstract
The maturation of inhibitory neurons is crucial for regulating plasticity in developing brains. Previous work has suggested that the Hurst exponent, the measure of autocorrelation in time series, reflects inhibition, but empirical data supporting this link are sparse. Here, we demonstrate significant spatial correlations between the Hurst exponent and ex vivo parvalbumin (PV) inhibitory mRNA expression in human children and adults, as well as between the Hurst exponent and PV-positive cell counts in mice, across both sexes. We further identified developmental plateaus in inhibition, as indicated by both PV inhibitory mRNA expression and the Hurst exponent, occurring prior to adolescence in humans and rats. In sum, this work suggests that the Hurst exponent can be used to study the development of inhibition in vivo and to understand inhibitory development across species.
Significance Statement
Understanding the spatial and temporal progression of developmental plasticity is crucial for identifying periods of vulnerability and opportunity in the human brain. This study supports the Hurst exponent as a promising in vivo marker for inhibition, showing significant spatial correlations with parvalbumin RNA expression and cell counts across species. Findings indicate a plateau in inhibition by late childhood in humans and prior to puberty in rats, suggesting the Hurst exponent's translational potential for exploring inhibitory development mechanisms in both humans and animals.
Introduction
Children's brains are highly plastic, enabling rapid learning and skill refinement. Converging evidence from human and animal neuroscience points to the importance of inhibitory cells for regulating plasticity (Hensch, 2005). These cells adjust timing and strength of synaptic inputs, critical for spike timing-dependent plasticity (Hensch, 2005). However, inhibition in the human brain has been hard to measure in vivo. As a result, the precise timing of changes in inhibition and its variation across the human cortex is unknown. Elucidating inhibitory development is key to understanding healthy brain maturation and identifying periods of vulnerability.
GABAergic neurons are classified as PV+ (∼40%), SST+ (∼30%), and VIP+ (∼10–15%; Xu et al., 2010; Almási et al., 2019). PV+ cells play a key role in plasticity by forming synapses on perisomatic regions of excitatory pyramidal cells, suppressing their spiking, resulting in shorter activation timescales to process fast-changing sensorimotor inputs (Wang et al., 2004; Y. Kim et al., 2017). Manipulations of PV+ cells lead to alterations in developmental plasticity (Fagiolini and Hensch, 2000; Di Cristo et al., 2007; van Versendaal and Levelt, 2016; Kalish et al., 2020). For example, accelerating the maturation of PV+ cells in the mouse visual cortex leads to an earlier onset and closure of the critical period for ocular dominance (Di Cristo et al., 2007). In contrast, destroying perineuronal nets around PV+ cells reactivates visual cortical plasticity in adult rats (Pizzorusso et al., 2002).
Postmortem histological studies have shown protracted PV+ maturation across species. PV levels plateau around adolescence in rats (Caballero et al., 2014; Gildawie et al., 2020), by 3–5 months in the macaque cortex (Condé et al., 1996) and around ages 5–10 years in the human prefrontal cortex (Grateron et al., 2003; Fung et al., 2010). However, pediatric postmortem data are limited, clinically biased and typically cover few brain areas.
Understanding human inhibitory development requires appropriate noninvasive in vivo tools. Magnetic resonance spectroscopy (MRS) can estimate GABA levels (Tognarelli et al., 2015), but findings on age-related changes have been inconsistent across prefrontal subregions and hemispheres (Perica et al., 2022). Furthermore, MRS can only sample one large voxel at a time, are motion-sensitive, and are not widely available in large developmental datasets (Kreis, 2004; Saleh et al., 2020; Porges et al., 2021).
An emerging alternative is the Hurst exponent, which quantifies long-range temporal correlations and scale-invariant dynamics in various time series (Linkenkaer-Hansen et al., 2001). Initially introduced to assess fractality, particularly the criticality of neuronal systems (Bak et al., 1987), the Hurst exponent has been shown to be reduced in numerous mental disorders and related traits, including attention-deficit/hyperactivity disorder (ADHD; Stier et al., 2023), schizophrenia (Sokunbi et al., 2014), depression (Wei et al., 2013), and substance abuse (Ide et al., 2016). Lombardi et al. demonstrated in vitro that a lack of inhibition significantly disrupts criticality, suggesting that inhibition is a fundamental feature of the critical state in neuronal networks (Lombardi et al., 2012). Building on these findings, computational models have demonstrated that enhancing inhibition leads to an increase in the Hurst exponent of simulated local field potential and BOLD signals (Gao et al., 2017; Lombardi et al., 2017; Trakoshis et al., 2020; Chini et al., 2022). Studies in rodents and nonhuman primates show that the Hurst exponent is sensitive to GABA changes induced by propofol or chemogenetic manipulation, supporting a causal link between GABA and the Hurst exponent (Gao et al., 2017; Trakoshis et al., 2020).
In the present study, we employed a cross-species, multimethod approach to characterize how the Hurst exponent changes across development. First, we investigated the spatial correlation between the Hurst exponent, calculated from resting-state fMRI data, and inhibitory marker mRNA expression in human children and adults. Next, we replicated this correlation in fMRI and cell density data from mice. We then traced the developmental trajectories of the Hurst exponent in humans, mice, and rats. Together, this work supports the use of the Hurst exponent to study the development of inhibition across the cortex.
Materials and Methods
Human data and analyses
RNA expression data
Publicly available RNAseq reads per kilobase per million (RPKM) data from the BrainSpan Atlas of the Developing Human Brain (https://www.brainspan.org/) were used to characterize patterns of inhibitory marker gene expression across development. RNAseq probes lacking Entrez IDs, which are unique numeric identifiers assigned to genes by the National Center for Biotechnology Information as part of the Entrez Gene database, were excluded (n = 30,984). To address duplicated probes, we retained only the probe with the highest mean expression, resulting in a dataset comprising a total of 21,377 genes, which were subsequently log2 transformed. We used the org.Hs.eg.db package in R to obtain the biological processes regulated by individual genes.
To validate the results, we also used gene expression data from the Allen Human Brain Atlas (Hawrylycz et al., 2012, 2015), which included postmortem samples of six donors (one female and five males) analyzed through microarray transcriptional profiling. Spatial maps of mRNA expression were provided in volumetric 2 mm isotropic MNI space, following enhanced nonlinear registration and whole-brain prediction using variogram modeling as implemented by Gryglewski et al. (2018).
fMRI data
The Institutional Review Board at the University of Pennsylvania approved all human fMRI research reported here. Adult data were from a previous study on individual differences in frontoparietal plasticity in humans (Boroshok et al., 2023). Prior to participation, all individuals provided informed, written consent. Participants were recruited through the University of Pennsylvania's study recruitment system, as well as via community and university advertisements. Inclusion criteria comprised proficiency in English, the absence of psychiatric or neurological disorders or learning disabilities, no recent or current illicit substance use, and no contraindications for MRI. A total of 61 participants completed MRI scans. Forty-six adult participants were included in the final sample (M, 21.39 years; SD, 1.91 years; 63% female; 39% male; 0% other or nonbinary; 24% Asian; 33% Black; 17% Hispanic/Latino; 4% Multiracial; and 19% White). The 77% of participants were undergraduates, and 18% were graduate students at the University of Pennsylvania. Exclusions included participants falling asleep during the scans (n = 3), recent illegal substance use not reported during screening but reported during participation (n = 1), inability to tolerate scanning (n = 1), and other technical issues (n = 10).
Child data were collected as part of a study on typical brain development. A subset of these data were previously published (n = 92; Tooley et al., 2022). All parents provided informed, written consent for their children's participation in the study. Verbal assent was obtained from children under 8 years old, while those 8 years old and older provided written assent. Recruitment efforts targeted Philadelphia and its surrounding areas, using advertisements on public transportation, partnerships with local schools, outreach programs, community family events, and social media ads. Children ranged in age from 4.11 to 10.59 years (M, 6.55; SD, 1.36; 47% male; 53% female; 0% other or nonbinary; 52.9% Black; 42.6% White; 12.5% Asian; and 0.5% others). The final child sample consisted of 128 participants. Exclusions included participants failing to complete the resting-state scan (e.g., due to falling asleep or opting to end the scan prematurely, n = 17) and instances where parents reported a diagnosis of ADHD or developmental delay during the visit, despite not disclosing it during screening (n = 4). Additionally, motion and quality criteria were applied to mitigate the impact of image quality on our analyses (see fMRI Image Quality and Exclusion Criteria section). These criteria were chosen to balance the need to maximize data inclusion in a young population (Leonard et al., 2017) with the necessity of minimizing the influence of low-quality data on connectivity metrics (Power et al., 2014).
For the replication analysis, we analyzed neuroimaging data from 630 healthy participants aged 5.58–21.92 years (54.8% female), obtained from the Human Connectome Project in Development (HCP-D; Release 2.0; Somerville et al., 2018).
fMRI image acquisition
Imaging was performed at the Center for Advanced Magnetic Resonance Imaging and Spectroscopy at the University of Pennsylvania. Scanning was conducted using a Siemens MAGNETOM Prisma 3 T MRI scanner with a Siemens 32-channel coil. Adult participants underwent a T1-weighted (T1w) multiecho (MEMPRAGE) structural scan [acquisition parameters, repetition time (TR), 2,530 ms; TI, 1,330 ms; TEs, 1.69/3.55/5.41/7.27 ms; BW, 650 Hz/px; 3× GRAPPA; flip angle, 7°; voxel size, 1 mm isotropic; matrix size, 256 × 256 × 176; field of view (FOV), 256 mm; total scan time, 4:38] using volumetric navigators and a 5 min resting-state run (TR, 2,000 ms; TEs, 30.20 ms; flip angle, 90°; resolution, 2 mm isotropic) with participants fixating on a cross throughout the scan. Resting-state scanning was continued until at least 5 min of data were acquired with framewise displacement <0.5 mm. Children first completed a mock scan to acclimate to typical MRI noises and practice remaining still. During the MRI session, a researcher remained in the scanner room to provide reassurance to the child. Children participated in the same structural and functional scans as adults. Head motion was monitored in real time using the Framewise Integrated Real-time MRI Monitor system, and 10 min of low-motion resting–state data (two runs with framewise displacement <1 mm) were collected when feasible. For participants with multiple usable resting-state runs, framewise displacement was averaged across runs, weighted by run length. All analyses were adjusted for average framewise displacement and the number of censored volumes.
In the HCP-D study, participants were instructed to remain still, stay awake, and blink normally while focusing on a fixation cross during the resting-state scan. For participants aged 8 years and older, 26 min of resting-state data were collected across four runs. For younger participants (ages 5–7 years), individual run durations were shortened, with a total scan time of 21 min.
fMRI image quality and exclusion criteria
The imaging data quality was evaluated through the fMRIPrep visual reports and MRIQC 0.14.2 software. Two reviewers manually evaluated all structural and functional images at each preprocessing stage for any image quality concerns. Functional images underwent visual inspection to ensure adequate whole-brain FOV coverage, the absence of signal blurring or artifacts, and correct alignment with the anatomical image.
For children's data, exclusions were made for participants with unusable structural images (n = 1), functional data artifacts (e.g., hair glitter, n = 1), incorrect scanner registration (n = 1), registration failure (n = 18), <100 s usable data after high-motion outlier volumes (n = 13), and average framewise displacement >1 mm (n = 6). For participants with multiple usable resting-state runs, framewise displacement was averaged across runs. To ensure that our results were not driven by motion, we conducted an additional analysis with a more stringent motion exclusion criterion. In this pipeline, we excluded participants who had average FD > 0.5 mm (n = 22), for a total of n = 108 participants. None of the adult data met the exclusion criteria of having an average framewise displacement >1 mm.
In the HCP-D study, participants with an average FD >1 mm were excluded (n = 51). For the sensitivity analysis, participants with an average FD >0.5 mm were excluded (n = 270).
fMRI image processing
The T1w image underwent intensity nonuniformity correction using N4BiasFieldCorrection (Tustison et al., 2010) from the Advanced Normalization Tools (ANTS) version 2.2.0 and used as T1w reference throughout the workflow. Following this, skull stripping was performed using the antsBrainExtraction.sh script (ANTS version 2.2.0), with the OASIS template as the target. The brain mask was refined with a custom variation of the method to reconcile ANTS-derived and FreeSurfer (Dale et al., 1999)-derived segmentations of the cortical gray matter of Mindboggle (RRID:SCR_002438). Spatial normalization to the ICBM 152 Nonlinear atlases version 2009c (RRID:SCR_008796) was performed through nonlinear registration with antsRegistration (ANTS, version 2.2.0, RRID:SCR_004757), using brain-extracted versions of both T1w volume and template. Brain tissue segmentation of CSF, white matter (WM), and gray matter was performed on the brain-extracted T1w using fast [fMRI of the Brain Software Library (FSL) version 5.0.9; RRID:SCR_002823].
For each resting-state BOLD run, a reference volume and its skull-stripped version were generated using a custom methodology of fMRIPrep version 20.2.0 (RRID:SCR_016216), which is based on Nipype version 1.1.7 (RRID:SCR_002502, nipy/nipype:1.1.7). Brain surfaces were reconstructed using the recon-all command before other processing, and reconstructed surfaces were used as input to fMRIprep. The BOLD reference was then coregistered to the T1w reference using bbregister (FreeSurfer), employing boundary-based registration (Greve and Fischl, 2009). Coregistration was configured with nine degrees of freedom to account for distortions remaining in the BOLD reference. Head-motion parameters with respect to the BOLD reference (transformation matrices and six corresponding rotation and translation parameters) were estimated before any spatiotemporal filtering using the mcflirt tool (FSL version 5.0.9), and slice-time correction was performed using 3dTshift from AFNI 20160207 (RRID:SCR_005927). The BOLD time series were resampled onto the MNIPediatricAsym standard space by applying a single, composite transform, generating a preprocessed BOLD run in MNI152NLin6Asym space.
Several confounding time series were calculated based on the preprocessed BOLD, including FD, DVARS (root-mean-square intensity difference from one volume to the next), and three region-wise global signals (CSF, WM, and the whole brain). FD and DVARS were calculated for each functional run, both using their implementations in Nipype (Power et al., 2014). The head-motion estimates calculated in the correction step were also placed within the corresponding confounds file.
All resamplings can be performed with a single interpolation step by composing all the pertinent transformations (i.e., head-motion transform matrices and coregistrations to anatomic and template spaces). Gridded (volumetric) resamplings were performed using antsApplyTransforms (ANTS), configured with Lanczos interpolation to minimize the smoothing effects of other kernels (Lanczos, 1964).
Further preprocessing was performed using a confound regression procedure that has been optimized to reduce the influence of participant motion; preprocessing was implemented in XCP-D 0.5.0 (Ciric et al., 2018), a multimodal tool kit that deploys processing instruments from frequently used software libraries, including FSL (Jenkinson et al., 2002) and AFNI (Cox, 1996). Further documentation is available at https://xcp-d.readthedocs.io/en/latest/ and https://github.com/PennLINC/xcp_d. Using XCP-D, the preprocessed functional time series on the fsLR cortical surface underwent nuisance regression, using a 36-parameter model that included three global signals (whole brain, CSF, and WM) and six motion parameters, as well as their derivatives, quadratic terms, and derivatives of quadratic terms. Linear regression, implemented in scikit-learn (0.24.2), was employed to regress the confounds. Motion censoring was applied, with outlier volumes exceeding FD = 0.3 mm flagged and removed from confound regression. An interpolated version of the BOLD data is created by filling in the high-motion outlier volumes with cubic spline interpolated data, as implemented in nilearn. Any outlier volumes at the beginning or end of the run are replaced with the closest nonoutlier volume's data, in order to avoid extrapolation by the interpolation function. If there was <100 s of data remaining after removing high-motion outlier volumes, those samples are excluded. No bandpass filtering was performed.
For HCP-D study, preprocessed structural MRI data were provided as part of the Lifespan HCP Release 2.0 through the NDA (https://nda.nih.gov/). The data were processed using HCP Pipelines (version 4.0.1) within the QuNex container environment (qunex.yale.edu). Detailed information of the minimal preprocessing pipelines can be found in Glasser et al. (2013). In brief, volumetric minimal preprocessing for structural MRI data includes gradient distortion correction, readout distortion removal, repetition averaging, and nonlinear registration to the MNI standard space. After completion of structural preprocessing, volumetric rs-fMRI data were minimally preprocessed, which comprised gradient distortion correction, head-motion correction, echoplanar imaging (EPI) distortion correction, denoising with ICA-FIX, native space and MNI standard space registration, and, lastly, global intensity normalization.
Hurst exponent
The Hurst exponent can be estimated using several methods. The power spectral density method leverages the relationship between the spectral density of the form 1/fβ and the Hurst exponent α, where β = 2α−1 (Ivanov et al., 2001). Alternatively, detrended fluctuation analysis (DFA) estimates the Hurst exponent by quantifying how fluctuations in a detrended, integrated time series scale with the size of the observation window, providing robustness against nonstationarities commonly present in neural data (Peng et al., 1994; Hu et al., 2001; Chen et al., 2005; Schmitt and Ivanov, 2007; Schmitt et al., 2009). Finally, fractional integration models treat the time series as fractionally integrated processes (FIP), where the Hurst exponent H is defined as H = d + 0.5, with d representing the fractional integration order (You et al., 2012; Trakoshis et al., 2020; Linke et al., 2024). In these models, the autocorrelation function (ACF) decays according to a power law, approximately as ACF(k) ∼ k2d−1, where k is the lag. This means that for d > 0, the correlations between distant observations decline slowly, reflecting long-range dependence or persistence in the signal. This approach fits parametric models to capture long-range dependencies more rigorously. Given the short and noisy nature of fMRI BOLD signals in developmental samples, we employed the fractional integration method in our analysis to ensure robust and stable estimation.
In this study, Hurst exponent was computed using the nonfractal MATLAB toolbox (https://github.com/wonsang/nonfractal; You et al., 2012; Trakoshis et al., 2020). The specific function utilized is bfn_mfin_ml.m function with the “filter” argument set to “haar” and the “ub” and “lb” arguments set to [1.5, 10] and [−0.5, 0], respectively. This toolbox uses a discrete wavelet transform and a model of the time series as FIP and is estimated using maximum likelihood estimation. The FIP model was employed instead of the fractional Gaussian noise model to eliminate the assumption of stationarity and the upper bound of H = 1. The FIP model provides greater flexibility and potentially increased sensitivity because it offers better estimation of between-subject variability, especially when estimates are close to or surpass H = 1. When H > 1, the time series is deemed nonstationary and displays long memory characteristics, such as being fractal.
Sensory–Association axis
We used the sensorimotor–association (S–A) axis derived by Sydnor et al. (2021). This map encompasses various cortical hierarchies, including functional connectivity gradients, evolutionary cortical expansion patterns, anatomical ratios, allometric scaling, brain metabolism measures, perfusion indices, gene expression patterns, primary modes of brain function, cytoarchitectural similarity gradients, and cortical thickness.
GAMs
To flexibly model linear and nonlinear relationships among the Hurst exponent and S–A axis rank, we employed generalized additive models (GAMs) using the mgcv package in R. GAMs, with the region-averaged Hurst exponent as the dependent variable, age as a smooth term, and gender, mean framewise displacement, and the number of censored as linear covariates. Models were fitted separately for each parcellated cortical region, using thin plate regression splines as the smooth term basis set and the restricted maximum likelihood approach for smoothing parameter selection. The smooth term for age generated a spline representing a region's developmental trajectory, with a maximum basis complexity (k) set to 3 to prevent overfitting.
For each regional GAM, we assessed the significance of the association between the Hurst exponent and age using an analysis of variance (ANOVA), comparing the full GAM model to a nested model without the age term. A significant result indicates that including a smooth term for age significantly reduced the residual deviance, as determined by the χ2 test statistic. For each regional GAM, we identified the specific age ranges where the Hurst exponent significantly changed using the gratia package in R. Age windows of significant change were determined by examining the first derivative of the age-smooth function (Δ Hurst exponent/Δ age) and assessing when the simultaneous 95% confidence interval of this derivative did not include 0 (two-sided). To quantify the overall magnitude and direction of the association between the Hurst exponent and age, referred to as a region's overall age effect, we calculated the partial R2 between the full GAM model and the reduced model for effect magnitude. We then signed the partial R2 based on the average first derivative of the smooth function for effect direction. We sorted 400 parcels based on their S–A rank and divided them into three groups: sensory, middle, and association regions. Subsequently, we performed the same GAM analysis for the averaged values across 133 parcels within each group.
Spin-based spatial permutation testing
To address distance-dependent spatial autocorrelation in cortical data, we evaluated the significance of Pearson's correlations between two whole-brain cortical feature maps (the Hurst exponent) using nonparametric spin tests (Alexander-Bloch et al., 2018). These tests, also known as spatially constrained rotation tests, compare the observed correlation to a null distribution obtained by spatially iterating (spinning) one of the feature maps. Specifically, the null distribution is generated by rotating spherical projections of one map while preserving its spatial covariance structure. The p value (pspin) is determined by comparing the empirical correlation to the distribution obtained from 10,000 spatial rotations. Spin tests were conducted using the spin permutation test algorithm available in the ENIGMA toolbox for Python (https://enigma-toolbox.readthedocs.io/en/latest/pages/08.spintest/index.html).
Mouse data and analyses
Inhibitory neuron cell density
For cell density analysis of PV+, SST+, and VIP+ cells, we used the dataset published by Y. Kim et al. (2017). Five male and five female transgenic mice were used per each cell type. Cell counting and distribution mapping were conducted using a platform developed by Y. Kim et al. (2017), which relies on automated imaging via serial two-photon tomography and data analysis through machine learning algorithms. Within this process, automatic cell counting within specific anatomical regions was achieved by tallying the number of registered signals within each anatomical label using custom-built codes. Cell counting was carried out in evenly spaced and partially overlapping 3D sphere voxels (with a diameter of 100 and 20 μm apart) to accurately digitize and visualize the distribution of detected signals throughout the entire brain in an unbiased manner.
fMRI data
For fMRI analysis, a publicly available fMRI dataset (doi:10.18112/openneuro.ds001653.v1.0.2) was used. All procedures performed in studies involving animals were in accordance with the ethical standards of the Institutional Animal Care and Use Committee (A*STAR Biological Resource Centre, Singapore, IACUC #161134/171203). For the analysis of correlation between cell density and the Hurst exponent, 10 female C57BL/6 mice (Janvier Labs), aged 3 months were included. A total of 57 runs were conducted, with each mouse undergoing six runs over two sessions. The mice were housed in standard mouse cages under a 12 h light/dark cycle, with access to food and water ad libitum. Anesthesia was induced with 4% isoflurane. Subsequently, the animals were placed on an MRI-compatible cradle, while isoflurane was reduced and kept to 1%. Care was taken to maintain the animal temperature at 37°C. Another 10 C57BL/6 female mice (Janvier Labs), aged 3 months, were included for testing the correlation under Med/Iso (medetomidine 0.05 mg/kg bolus/0.1 mg/kg/h infusion, i.v., +isoflurane 0.5%) anesthesia protocol. A total of 57 runs were conducted, with each mouse undergoing six runs over two sessions. For those mice, anesthesia was induced with 4% isoflurane. Subsequently, a bolus with a mixture of medetomidine (Dormitor, Elanco) and pancuronium bromide (muscle relaxant, Sigma-Aldrich) was administered subcutaneously (0.05 mg/kg), followed by a maintenance infusion (0.1 mg/kg/h), while isoflurane was simultaneously reduced and kept to 0.5%.
To analyze the developmental trajectory of the Hurst exponent, we included a total of 72 wild-type (comprising 40 females and 32 males) within the age range of 1–8 months, distributed across four distinct age groups: 1 month (n = 21, with 9 females and 12 males), 2 months (n = 16, with 6 females and 10 males), 5 months (n = 17, with 13 females and 4 males), and 8 months (n = 18, with 12 females and 6 males). The mice were housed in standard mouse cages under a 12 h light/dark cycle, with access to food and water ad libitum. Anesthesia was induced using isoflurane (Abbott) in a 1:4 oxygen to air mixture, with a concentration of 3% for induction and 1.4% during measurements, administered via a facemask. During the experiments, the mice were positioned on a mouse support apparatus with ear bars used to minimize motion. Body temperature was continuously monitored using a rectal thermometer and maintained at 36.5 ± 0.5°C through an adjustable warm water bath integrated into the support. The entire procedure, including mouse preparation, typically lasted for 45 min.
fMRI image acquisition
Data were acquired on an 11.75 T scanner (Bruker BioSpin MRI) equipped with a BGA-S gradient system, a 72 mm linear volume resonator coil for transmission, and a cryoprobe 2 × 2 phased-array surface coil. Images were acquired using the Paravision 6.0.1 software. An anatomical reference scan was acquired using a spin-echo Turbo-RARE sequence: FOV, 17 × 9 mm2; FOV saturation slice masking nonbrain regions; number of slices, 28; slice thickness, 0.35; slice gap, 0.05 mm; matrix dimension (MD), 200 × 100; TR, 2,742 ms; echo time (TE), 30 ms; RARE factor, 8; and number of averages, 2. fMRI was acquired using a gradient-echo EPI (GE-EPI) sequence with the same geometry as the anatomical scan: MD, 90 × 60; in-plane resolution, 0.19 × 0.15 mm2; TR, 1,000 ms; TE, 11.7 ms; flip angle, 50°; volumes, 180; and bandwidth, 119,047 Hz. Field inhomogeneity was corrected using the MAPSHIM protocol.
fMRI image processing
A study-specific template was created from the individual anatomical images using the buildtemplateparallel.sh script from the ANTs (20150828 builds). The study template was subsequently registered to the reference template of the Allen Mouse Common Coordinate Framework version 3 (CCF v3; Oh et al., 2014). Individual anatomical images were corrected for b1 field inhomogeneity (N4BiasFieldCorrection), denoised (DenoiseImage), intensity thresholded (ImageMath), brain masked (antsBrainExtraction.sh), and registered to the study template, resampled to a 0.2 × 0.2 × 0.2 mm3 resolution, using a SyN diffeomorphic image registration protocol (antsRegistration). The functional EPI time series were despiked (3dDespike) and motion-corrected (3dvolreg) using the Analysis of Functional NeuroImages software (AFNI_16.1.15). A temporal average was estimated using fslmaths (FMRIB Software Library 5.0.9), corrected for the b1 field (N4BiasFieldCorrection), brain masked (antsBrainExtraction.sh), and registered to the anatomical image (antsRegistration). The fMRI time series were bandpass filtered (3dBandpass) to a 0.01–0.25 Hz band. An independent component analysis (MELODIC) was performed and used as a basis to train a FIX classifier (Salimi-Khorshidi et al., 2014) in 15 randomly selected runs. The classifier was applied to identify nuisance-associated components that were regressed from the fMRI time series. The EPI to anatomical, anatomical to study template, and study template to Allen Mouse CCF v3 registrations were combined into one set of forward and backward transforms (ComposeMultiTransform). The combined transforms were applied to the denoised fMRI data. Signal time series for each region of interest (ROI) in the Allen Mouse Brain atlas (Goldowitz, 2010) were extracted.
Hurst exponent
We used the same metrics as those applied to humans to calculate the Hurst exponent for mice.
Rat data and analyses
fMRI data
For the fMRI analysis, we utilized the Standard Rat dataset as described by Grandjean et al. (2023). All fMRI acquisitions were conducted following approval from the relevant local and national ethics authorities. To investigate the developmental trajectory of the Hurst exponent, we included a total of 145 rats, comprising 54 females and 91 males, aged between 1 and 4 months. These rats were distributed across four distinct age groups: 1.0 month (n = 10, males only), 1.5 months (n = 10, males only), 2.0 months (n = 68, with 26 females and 42 males), 2.5 months (n = 20, with 5 females and 15 males), 2.9 months (n = 3, females only), 3.0 months (n = 14, with 11 females and 3 males), and 4.0 months (n = 20, with 9 females and 11 males).
fMRI image acquisition
Rats were anesthetized using 4% isoflurane and a subcutaneous bolus of 0.05 mg/kg medetomidine for induction, followed by maintenance anesthesia with 0.4% isoflurane and a continuous infusion of 0.1 mg/kg/h medetomidine subcutaneously. Imaging was conducted using a GE-EPI technique 40 min after anesthesia induction. The imaging parameters included a TR of 1,000 ms, TE/flip angle/bandwidth defined according to field strength, 1,000 repetitions, a matrix size of 64 × 64, a FOV of 25.6 × 25.6 mm2, and 18 interleaved axial slices of 1 mm thickness with a 0.1 mm gap. The full protocol can be found at https://github.com/grandjeanlab/StandardRat.
fMRI image processing
The scans were organized following the BIDS format. Preprocessing was conducted separately for each scan session using a reproducible containerized software environment for RABIES 0.3.5 (Singularity 3.7.3–1.el7, Sylabs). This involved several steps: autobox, N4 inhomogeneity correction, motion correction, rigid registration between functional and anatomical scans, nonlinear registration between anatomical scan and template, and resampling to a common space at 0.3 × 0.3 × 0.3 mm3 resolution. To address differences in the brain size, image contrast, and susceptibility distortions in rodent images, we developed a workflow for volumetric image registration and brain masking. Quality control was ensured through visual inspection of preprocessing outputs for all scans. Global signal regression was performed along with motion regression, spatial smoothing to 0.5 mm3, a high-pass filter at 0.01 Hz, and a low-pass filter at either 0.1 or 0.2 Hz. The ICA-AROMA method was adapted from humans to rats by utilizing dedicated rat cerebrospinal fluid and brain edge masks and by training the classifier parameters based on a set of rodent images. Visual inspection of the components and their classifications showed that <5% of the plausible signal components were mistakenly labeled as noise. Signal time series for each ROI in the SIGMA anatomical atlas (Barrière et al., 2019) were extracted using the NiftiLabelsMasker function in Nilearn 0.7.1.
Hurst exponent
We used the same metrics as those applied to humans to calculate the Hurst exponent for rats.
All the code is publicly available at https://github.com/monami-nishio/hurst_exponent_development/tree/main.
Results
Hurst exponent correlates with parvalbumin RNA expression in humans
To assess the validity of the Hurst exponent as an inhibitory marker, we investigated the correlation between the Hurst exponent and mRNA inhibitory neuronal markers from the BrainSpan Atlas of the Developing Human Brain (BrainSpan Atlas) and Allen Human Brain Atlas (Hawrylycz et al., 2012; Gryglewski et al., 2018). To calculate the Hurst exponent, first we parcellated the cortex into 400 equally sized regions using Schaefer 400 parcellation (Schaefer et al., 2018). Using fractional integration models (You et al., 2012; Trakoshis et al., 2020; Linke et al., 2024), we computed the ACF of the time series for each cortical parcel, estimated the fractional integration order (d), and the Hurst exponent calculated as H = d + 0.5 (Fig. 1A). The Hurst exponent across the cortex shows a U-shaped pattern along the S–A axis, with both the sensory and association ends exhibiting higher values than the middle regions, in both children (Fig. 1C,E; ages 4–11 years, n = 128, R2 = 0.037; p = 0.018) and adults (Fig. 1D,F; ages 18–25 years, n = 46, R2 = 0.050; p < 0.001). As expected, adults exhibited significantly higher Hurst exponents compared with children throughout the cortex, suggesting elevated levels of inhibition in adults relative to children (children vs adults, independent t test, t(172) = 4.301; p < 0.001). As mRNA inhibitory neuronal markers, we analyzed three markers of specific inhibitory neuron subtypes: parvalbumin (PV), somatostatin (SST), and vasoactive intestinal polypeptide (VIP; Fig. 1B). We also analyzed the expression level of glutamic acid decarboxylases (GAD) 1 and 2, which catalyze the production of GABA and are expressed across all types of inhibitory neurons (Fig. 1B).
The Hurst exponent across the cortex. A, Examples of fMRI time series with high and low Hurst exponents and ACF. ACF decays according to a power law, approximately as ACF(k)∼ k2d−1, where k is the lag. H is defined as H = d + 0.5, with d representing the fractional integration order. B, Illustration of inhibitory neuron types, including PV, SST, VIP, as well as proteins involved in GABA production, glutamic acid decarboxylase (GAD 1 and 2). C, D, Group-averaged Hurst exponent maps for children (ages 4–11 years, n = 128; C) and adults (ages 18–25 years, n = 46; D). E, F, Mapping of the Hurst exponent along the S–A axis for children (E) and adults (F). The GAM-predicted Hurst exponent value at each parcel along the S–A axis is displayed with a 95% confidence interval (R2 = 0.037; p = 0.018 for children; R2 = 0.050; p < 0.001 for adults).
For the BrainSpan Atlas, we included mRNA data from all 11 Brodmann's areas (Brodmann, 1909) sampled at ages 3, 8, 11, and 13 years for children (n = 6) and ages 18, 19, 21, and 23 years for adults (n = 4; Fig. 2A). As there were only two regions sampled for a 4-year-old, we excluded these data. No data were available for other ages. We aligned the Brodmann's areas with the Schaefer 400 parcellation to extract Hurst exponents and to assign the regions to an S–A axis rank.
The Hurst exponent correlates with PV mRNA expression across the cortex. A, Brodmann's areas from which inhibitory marker expression was analyzed. Brain regions are arranged along a spectrum from sensorimotor (yellow) to association (purple). Vis, visual cortex (BA17); Aud, auditory cortex (BA41); Somat, somatosensory cortex (BA1–3); Mot, motor cortex (BA4); STG, superior temporal gyrus (BA22); IPL, inferior parietal lobe (BA39, 40); OFC, orbitofrontal cortex (BA11–14); ACC, anterior cingulate cortex (BA24, 25); IT, inferior temporal gyrus (BA20); vlPFC, ventrolateral prefrontal cortex (BA44, 45); dlPFC, dorsolateral prefrontal cortex (BA8, 9, 46). B, C, Correlation between PV expression (Log2 RPKM mapped reads) and the Hurst exponent in children (B; r = 0.641; pspin = 0.034; rpartial = 0.827; ppartial = 0.003) and adults (C; r = 0.695; pspin = 0.018; rpartial = 0.778; ppartial = 0.008). PV expression was extracted from Brodmann's areas, and Hurst exponents were extracted from the Schaefer-400 parcels with >50% overlap with each Brodmann's area. rpartial and ppartial represent the results of the partial correlation test, accounting for the S–A hierarchy as a covariate. D, Geodesic distance (mm) from primary sensory centers, which were defined as the geometric centers of the primary somatosensory, visual, and auditory areas (Huntenburg et al., 2021). Brain regions are arranged along a spectrum from closer (green) to farther (pink) from the primary sensory centers. Anatomical segmentations are based on the Allen Mouse Brain atlas. E, Correlation between PV+ cell density and Hurst exponent (r = 0.790; pspin = 0.001). Sensory regions (yellow) include VISp, primary visual area; AUDp, primary auditory area; SSp-n, primary somatosensory area, nose domain; SSp-m, primary somatosensory area, mouth domain; SSp-bfd, primary somatosensory area, barrel field; and SSs, supplementary somatosensory area. Motor regions (orange) include Mop, primary somatomotor area, and MOs, secondary somatomotor area. Association regions (purple) include TEa, temporal association areas; ECT, ectorhinal area; PERI, perirhinal area; AI, agranular insular area; RSP, retrosplenial area; ACA, anterior cingulate area; and PL, prelimbic area.
Using the BrainSpan Atlas, the distribution of the Hurst exponent significantly correlated with the PV mRNA expression both in children (Fig. 2B; r = 0.641; pspin = 0.034) and in adults (Fig. 2C; r = 0.695; pspin = 0.018). Sensory and motor regions tended to have higher PV expression, as previously established (Anderson et al., 2020), as well as higher Hurst exponents. In contrast, the Hurst exponent was negatively correlated with SST expression and showed no correlation with VIP, GAD1, or GAD2 expression in neither children or adults (Fig. S1).
Given that PV expression is notably higher in sensory regions compared to association regions, one possibility is that the correlation between the Hurst exponent and PV expression stems solely from their alignment along the S–A axis. If this is the case, we would expect similarly strong correlations between other genes aligned with the S–A axis and the Hurst exponent. We found 655 genes significantly correlated with the S–A axis in both children and adults (p < 0.05). However, of these 655 genes, only four were significantly correlated with the Hurst exponent (PVALB, KCNS1, JDP2, ITIH5). Among these, PV had a numerically higher R2 for both children and adults’ Hurst exponent compared with the other three genes. We also performed a partial correlation test and found that even when controlling for S–A axis rank, PV expression and the Hurst exponent still showed a significant correlation (children, rpartial = 0.827; ppartial = 0.003; adults, rpartial = 0.778; ppartial = 0.008). This suggests that their correlation cannot be solely attributed to their alignment with the S–A hierarchy.
Another possibility is that although PV significantly correlates with the Hurst exponent, other genes might correlate even more strongly and better explain the spatial distribution of the Hurst exponent. Out of the 21,377 unique genes in our dataset, only 42 genes (<0.2%) displayed significant correlations (p < 0.05) with both child and adult Hurst exponents. Among these 42 genes, only eight other genes (FNDC4, PIM3, F10, SNX21, UBE2QL1, SEZ6L2, RXRB, IL32) had a numerically higher R2 for both children and adults than PV. These eight genes encode proteins known to play key roles in metabolism.
Since PV expression is activity-dependent and PV interneurons have exceptionally high metabolic demands (Donato et al., 2013), the fact that these eight metabolism-related genes are more correlated with the Hurst exponent than PV raises an important question: Does the Hurst exponent truly track PV-positive cell density or is it instead a proxy for overall regional activity and metabolic load? To test the “metabolism” hypothesis directly, we compared the adult Hurst exponent map with established cortical metabolism maps. Specifically, we examined cerebral blood flow (CBF), cerebral blood volume (CBV), oxygen consumption (CMRO₂), and glucose utilization (CMRglc) derived from the PET study (Vaishnavi et al., 2010). While CBF showed a slight positive correlation with the Hurst exponent (rpartial = 0.16; ppartial = 0.001), CBV showed a significant negative correlation (rpartial = −0.1; ppartial = 0.045), and neither CMRO₂ nor CMRglc was significantly correlated (CMRO₂, rpartial = −0.01; ppartial = 0.816; CMRglc, rpartial = 0.05; ppartial = 0.367). Taken together, these additional analyses do not support an alternative hypothesis that the Hurst exponent reflects global activity or metabolic levels across cortical regions.
For replication, we used the Allen Human Brain Atlas to test the correlation between PV mRNA and the Hurst exponent. We parcellated the whole-brain mRNA expression maps (Gryglewski et al., 2018) using the Schaefer 400 parcellation (Schaefer et al., 2018). With higher resolution, we successfully replicated the significant correlation between PV mRNA expression and the Hurst exponent (Fig. S2F). Additionally, with this enhanced resolution, we identified significant positive correlations between GAD1 and GAD2 mRNA expression and the Hurst exponent (Fig. S2I,J). All correlations remained significant even after accounting for the S–A axis using a partial correlation test.
Hurst exponent correlates with PV cell density in mice
Next, we assessed whether the correlation between the Hurst exponent and PV+ cells extended to mice. In lieu of mRNA expression data, we used cell density (cells per cubic millimeter) for PV+, SST+, and VIP+ cells in mice (Y. Kim et al., 2017). For the Hurst exponent calculation, we used a publicly available dataset of mouse resting-state fMRI (doi:10.18112/openneuro.ds001653.v1.0.2). The Allen Mouse Brain Atlas was used for anatomical parcellation of both cell density and fMRI whole brain map. We arranged brain regions of the Allen Mouse Brain Atlas based on the geodesic distance from primary sensory centers, which were defined as the geometric centers of the primary somatosensory, visual, and auditory areas (Fig. 2D, reproduced from Huntenburg et al., 2021). Geodesic distance measures the shortest path between two points on a surface, following the surface's contours (Huntenburg et al., 2021). Distances were calculated using the surface mesh representation of the Allen Mouse CCF v3 (Sanzleon et al., 2013).
As in humans, we found relatively higher PV+ cell density and higher Hurst exponents in sensory and motor regions compared with association regions (Fig. 2E). The Hurst exponent was strongly positively correlated with the density of PV+ cells (Fig. 2E; r = 0.790; pspin = 0.001). The Hurst exponent was negatively correlated with SST+ and did not correlate with VIP+ cells (Fig. S1I,J). We did find that the anesthesia protocol significantly influences the correlations between the Hurst exponent and cell counts. We observed a significant positive correlation between PV+ cell density and the Hurst exponent in mice anesthetized with 1% isoflurane, but this correlation was eliminated under the Med/Iso (medetomidine 0.05 mg/kg bolus/0.1 mg/kg/h infusion, i.v., +isoflurane 0.5%) anesthesia protocol (r = −0.161; pspin = 0.278).
Plateaus in inhibitory development across humans, mice, and rats
We first analyzed the relationship between mRNA inhibitory markers and age during childhood in humans across 11 cortical regions, using data from the BrainSpan Atlas (Fig. 3A). We observed increases in mRNA expression which plateaued by the age of 8 for all inhibitory markers (Fig. 3B–F; all ages greater than age 3; ANOVA F(3,55) > 23.531; p < 0.001). Additionally, the expression of PV, SST, and VIP at age 13 was significantly lower compared with age 8 (post hoc Tukey's HSD p's < 0.05), whereas no significant differences were found for GAD1 and GAD2. Associations with age did not differ across regions of the cortex (p's > 0.05).
Developmental plateaus of inhibition. A, Brodmann's areas (BA) from which inhibitory marker expression was analyzed. Brain regions are arranged along a spectrum from sensorimotor (yellow) to association (purple). B–F, Expression of PV (B, ANOVA F(3,55) = 23.531; p < 0.001), SST (C, ANOVA F(3,55) = 99.770; p < 0.001), VIP (D, ANOVA F(3,55) = 31.053; p < 0.001), GAD1 (E, ANOVA F(3,55) = 301.203; p < 0.001), and GAD2 (F, ANOVA F(3,55) = 197.496; p < 0.001) across development from age 3–13 years. Line colors correspond to colors in A. ***p < 0.001, post hoc Tukey HSD test comparing each older age group with 3-year-olds. Error bars indicate ±SD across donors. G, The 400 parcels are segmented along the S–A axis into three bins; sensory (yellow), middle (orange), and association (purple). Each bin contains 133 parcels. H, GAM-predicted developmental trajectories of the Hurst exponent in child dataset. Regional trajectories for each bin illustrate the GAM-predicted Hurst exponent values at each age, accompanied by a 95% confidence interval (sensory, partial R2 = 0.016; pFDR < 0.001; middle, partial R2 = 0.016; pFDR < 0.001; association, partial R2 = 0.012; pFDR < 0.001). A dashed line marks the age at which the first derivative of the age-smooth function (Δ Hurst exponent/Δ age) becomes insignificant. The age distribution of participants is depicted as a histogram along the x-axis. I, GAM-predicted developmental trajectories of the Hurst exponent in HCP-D dataset. Regional trajectories for each bin illustrate the GAM-predicted Hurst exponent values at each age, accompanied by a 95% confidence interval (sensory, partial R2 = 0.010; pFDR < 0.001; middle, partial R2 = 0.013; pFDR < 0.001; association, partial R2 = 0.012; pFDR < 0.001). J, The Hurst exponent across development of mice (sensory, ANOVA F(3,107) = 1.237; p = 0.298; association, ANOVA F(3,107) = 0.961; p = 0.431). Sensorimotor regions are shown in yellow, while association regions are shown in purple. K, The Hurst exponent across development of rats (sensory, ANOVA F = 4.806; p = 0.005; association, ANOVA F(6,138) = 2.437; p = 0.068). **p < 0.01, the post hoc Tukey HSD test comparing each older age group with 1-month-olds. Error bars indicate ±SD across subjects. Sensorimotor regions are shown in yellow, while association regions are shown in purple.
Next, we investigated whether a plateau in inhibition was also detectable in the Hurst exponent in children. We fitted region-specific GAMs with a smooth term for age. Each GAM estimated a smooth function (the model age fit) that described the relationship between the Hurst exponent and age, thereby modeling a region's developmental trajectory. Gender, mean framewise displacement, and the number of censored volumes were included as covariates; slopes did not differ by gender. We calculated the partial R2 change between the full GAM model and the reduced model (no age effects) and signed the partial R2 by the sign of the average first derivative of the smooth age function (effect direction). The 400 parcels were divided along the S–A axis into three bins, each containing 133 parcels (Fig. 3G). Across sensory, middle, and association regions, the Hurst exponent plateaued ∼8 years old [Fig. 3H, sensory, partial R2 = 0.016; pFDR < 0.001; middle, partial R2 = 0.016; pFDR < 0.001; association, partial R2 = 0.012; pFDR < 0.001; a dashed line indicates the age at which the first derivative of the age-smooth function (Δ Hurst exponent/Δ age) becomes insignificant]. Age effects did not differ across S–A bins.
For replication, we used the HCP-D (Somerville et al., 2018) dataset to track the developmental trajectory of the Hurst exponent across a wider age range. For sensory regions, the plateau timing was around 8 years old, replicating the previous result (Fig. 3I, sensory, partial R2 = 0.010; pFDR < 0.001). In contrast, the middle and association regions exhibited slightly delayed plateaus, occurring ∼10 years old (Fig. 3I, middle, partial R2 = 0.013; pFDR < 0.001; association, partial R2 = 0.012; pFDR < 0.001). After reaching a plateau in late childhood, the Hurst exponent declined across all sensory, middle, and association regions during adolescence and early adulthood.
To understand the developmental trajectory of the Hurst exponent in mice, we used the fMRI data obtained at 1, 2, 5, and 8 months old, roughly corresponding to adolescence and young adulthood (Dutta and Sengupta, 2016). We divided the cortex into sensorimotor regions (eight sensory regions, two motor regions) and association regions (seven regions). In aggregate, sensorimotor and association regions did not show significant age effects in this age range (Fig. 3J, sensorimotor, ANOVA F(3,107) = 1.237; p = 0.298; association, ANOVA F(3,107) = 0.961; p = 0.431).
Since the developmental stage most similar to human childhood occurs before 1 month in mice, and mice younger than 1 month are difficult to image with fMRI, we conducted an additional developmental analysis using rat fMRI. We used a public rat fMRI dataset with data available for 1.0, 1.5, 2.0, 2.5, 2.9, 3.0, and 4.0 months (Grandjean et al., 2023), an age range that roughly corresponds to peripubertal (1.0, 1.5 months), adolescence (2.0 months), and adulthood (2.5–4.0 months; Ghasemi et al., 2021). The SIGMA anatomical atlas (Barrière et al., 2019) was employed for whole-brain parcellation. Within this age range, sensory regions showed a significant age effect, while association regions exhibited a trending age effect (Fig. 3K, sensorimotor, ANOVA F(6,138) = 4.806; p = 0.005; association; ANOVA F(6,138) = 2.437; p = 0.068). The post hoc Tukey HSD test revealed that the Hurst exponent in sensorimotor regions was significantly higher at 1.5 and 2.0 months compared with 1.0 month (1.5 month; Q(138) = −0.032; p = 0.019; 2.0 months; Q(138) = −0.028; p = 0.003). No significant differences were observed between 1.0 month and ages beyond 2.5 months (2.5 months; Q(138) = −0.020; p = 0.215; 2.9 month; Q(138) = −0.028; p = 0.427; 3.0 month; Q(138) = −0.023; p = 0.149; 4.0 month; Q(138) = 0.016; p = 0.494).
Sensitivity analyses
We conducted several sensitivity analyses to confirm the robustness of our results. First, we conducted our main analyses with more stringent motion exclusion criteria (mean framewise displacement <0.5 mm). No adults were excluded under this stricter criterion. However, 22 participants were excluded from the analyses of the children's data, and 270 participants were excluded from the analyses of the HCP-D data. Even under the stringent motion exclusion criteria, the relationship between PV expression and the Hurst exponent was significant in children (rpartial = 0.801; pspin = 0.005). Developmental trajectories followed a similar shape (ps < 0.001 for sensory, middle, and association regions), but the timing of the plateau appeared slightly different in both child data (Fig. S3A) and HCP-D data (Fig. S3B).
Second, we replicated our key findings using an alternative method for estimating the Hurst exponent—DFA. The developmental trajectory of the Hurst exponent was consistently replicated across both child data and HCP-D data (Fig. S3C,D; ps < 0.001 for sensory, middle, and association regions). Additionally, this approach revealed even stronger correlations between PV mRNA expression and the Hurst exponent in human children (Fig. S3E; rpartial = 0.864; pspin = 0.001), as well as between mouse PV cell density and the Hurst exponent (Fig. S3F; r = 0.811; pspin < 0.001).
Discussion
We employed a cross-species, multimodal approach to investigate the early development of inhibitory circuits. We found significant spatial correlations between the Hurst exponent and the expression of PV mRNA across cortex in humans and in mice. We also found that PV expression and the Hurst exponent show plateaus by late childhood across sensorimotor and association areas. These results support the use of the Hurst exponent as a noninvasive method to investigate the development of inhibition and, by extension, plasticity in humans.
In this study, we discovered a correlation between the Hurst exponent and the expression of PV mRNA in the human cortex. We also found that the Hurst exponent correlates with the density of PV+ cells across the cortex in mice. In both species, the Hurst exponent was more correlated with PV cells than other inhibitory subtypes such as SST+ and VIP+ cells. This variation between inhibitory subtypes might be explained by the different functions of each subtype (Kepecs and Fishell, 2014). PV+ cells predominantly form synapses on perisomatic regions, exerting significant influence over the output of pyramidal neurons (Kepecs and Fishell, 2014). In contrast, SST+ cells primarily target dendrites to regulate inputs, while VIP+ cells synapse with both PV+ and SST+ cells, thereby not directly regulating either the input or output of pyramidal cells (Chaudhuri et al., 2015; C. K. Kim et al., 2017).
It is important to note that the expression of PV mRNA is activity-dependent (Donato et al., 2013). While we show that the Hurst exponent significantly correlates with PV cell counts in mice, we also found that brain state—specifically the anesthesia protocol—significantly influences the relationship between inhibitory neuron density and the Hurst exponent. Similarly, in humans, the Hurst exponent is known to vary with cognitive state and has been shown to correlate with task performance (Suckling et al., 2008; Kardan et al., 2023). Given these state-dependent effects on the Hurst exponent, further research is needed to clarify precisely what the Hurst exponent reflects in relation to PV neurons—whether it indexes cell count, neural activity, or other related properties.
The Hurst exponent is also widely recognized as a measure of long-range, power-law (scale-invariant) autocorrelation. A previous study provided cross-species evidence of such scale-invariant temporal patterns in brain dynamics during sleep–wake transitions, supporting the idea that criticality may be a conserved organizational principle across mammalian brains (Lo et al., 2004). One possible explanation for why PV-mediated inhibition is associated with a higher degree of long-range autocorrelations lies in the role of PV+ interneurons in regulating the balance and timing of cortical excitation. These fast-spiking interneurons enable precise temporal control of neuronal firing (Wang et al., 2004; Y. Kim et al., 2017). Enhanced PV+ inhibition may suppress high-frequency, uncorrelated (noise-like) activity, thereby promoting more coherent, low-frequency fluctuations that are reflected in higher Hurst exponents. This shift toward more structured activity could indicate a network operating closer to criticality—where brain dynamics are neither random nor overly synchronized—supporting both flexibility and stability in information processing (Bak et al., 1987). In this sense, increased PV-mediated inhibition may help support persistent activity and long-range dependencies, consistent with scale-invariant dynamics. This interpretation aligns with prior work highlighting the importance of inhibitory balance in sustaining critical brain states (Lombardi et al., 2012, 2017).
Given that PV+ cells are the critical regulators of plasticity in the developing brain (Fagiolini and Hensch, 2000; Di Cristo et al., 2007; van Versendaal and Levelt, 2016; Kalish et al., 2020), this finding highlights the potential of the Hurst exponent as a promising in vivo marker for measuring children's brain plasticity. We observed a plateau in inhibition by late childhood in both the mRNA and fMRI human data, as well as in rat fMRI data. This result is consistent with a previous study examining postmortem human samples of the dorsolateral prefrontal cortex, which noted an increase in PV gene expression from neonates to late childhood (Fung et al., 2010). It is also consistent with MRS findings showing an increase in GABA levels during childhood, followed by a plateau from adolescence through midlife, integrating results from multiple studies (Porges et al., 2021). In the human primary visual cortex, another metric of inhibition, the balance of GABAAα1:GABAAα2 receptor subunits reaches the adult level at the age of 4–5 years (Pinto et al., 2010). The expression level of gephyrin, a GABAA receptor anchoring protein, also shows a peak in expression during childhood (5–11 years old; Pinto et al., 2010). These results are further supported by behavioral studies showing that skills such as reasoning, executive function, emotion regulation, and social cognition show significant improvement until ∼8–10 years old, after which the rate of change decreases (Akshoomoff et al., 2014; Wellman, 2014; Fortenbaugh et al., 2015; Whitaker et al., 2018). Our finding that inhibition plateaus by late childhood suggests that this stabilization of inhibition may contribute to the slowdown in skill acquisition.
Interestingly, we also observed a subsequent decline following the late childhood plateau in human PV, SST, and VIP mRNA, as well as in human and rat Hurst exponents. This finding aligns with a previous MRS study that reported a decrease in GABA/Cr levels in the human anterior cingulate cortex during adolescence (Perica et al., 2022). This reduction in inhibition could be explained as an adjustment after an overshoot of inhibition, similar to processes like synaptogenesis and pruning, or as an adjustment in inhibitory activity to maintain the excitatory/inhibitory (E/I) balance following the reduction of excitatory activity during synaptic pruning in adolescence. Future research is needed to explore the underlying biological mechanisms driving these inhibitory changes into adolescence and early adulthood.
Our findings on maturational differentiation along the S–A axis were mixed. In the mRNA and our child dataset, we found no differences in the timing of the Hurst exponent plateau across the S–A axis, with all regions plateauing ∼8 years old. However, in the HCP-D dataset, middle and association regions showed a delayed plateau, occurring ∼10 years old, compared with sensory regions, which plateaued ∼8 years old. The mixed results might stem from the age range studied, as including younger or older individuals could alter the model's fit (Sydnor et al., 2021, 2023; Larsen et al., 2023). Given the smaller sample size in our child dataset after age 8, the HCP-D estimates for middle and association regions may be more reliable. For a comprehensive understanding of inhibitory development throughout the lifespan, datasets covering all age ranges from early infancy to adulthood, with balanced sample sizes across age groups, are essential.
Several potential limitations should be noted. First, the mRNA dataset for humans only included 1–2 children per age and did not include children between the ages of 5 and 7 years. Additionally, the number of brain regions included was limited. Future work with longitudinal data in a broader age range and including denser sampling of brain regions will be necessary to refine the temporal and spatial sequences of the relationships we report, as well as to better evaluate nonlinearities. Second, by carefully excluding data with motion artifacts, we may have limited the generalizability of our findings. Most young children move in the scanner, so it is essential to develop more motion-resilient sequences to allow investigators to acquire data from a more representative sample of young kids. Third, the age range of rodent datasets is limited because mice under 30 d old, corresponding to childhood, cannot be scanned due to their small brain size.
In summary, our study builds support for the validity of the Hurst exponent as an inhibitory marker across species. Understanding the developmental processes of inhibition in the human brain is crucial for characterizing both typical and atypical brain development. Inhibitory activity is reduced in autism (Trakoshis et al., 2020), schizophrenia (Liu et al., 2021), and depression (Gerner and Hare, 1981). Additionally, early life adversity, including maltreatment (Hanson and Nacewicz, 2021) and prenatal exposure to substances such as alcohol (Shenoda, 2017) and cocaine (McCarthy et al., 2011) can impact inhibitory development. By employing the Hurst exponent as an in vivo marker of inhibition in the developing human brain, we can investigate mechanisms of inhibitory development and guide interventions for risk factors associated with atypical inhibitory development.
Footnotes
First, we thank all of the individuals who participated in this research. We also thank Maayan Ziv, Steven Lopez, Cassandra Garcia, Priya Deliwala, Jasmine Forde, Katrina Simon, Sophie Sharp, Yoojin Hahn, Stephanie Bugden, Jamie Bogert, Alexis Broussard, Ava Cruz, Samantha Ferleger, Destiny Frazier, Jessica George, Abigail Katz, Sun Min Kim, Hunter Liu, Dominique Martinez, Ortal Nakash, Emily Orengo, Christina Recto, Leah Sorcher, Alexis Uria, Adam Pines, Quentin Wedderburn, Connor Kendzora, Anna Meaney, Leah Sorcher, Aparna Ramanujam, Daniel Southwick, Alexis Broussard, and Andrea Gomez for their help with data acquisition and cleaning. HCP-D data were obtained from the National Institute of Mental Health (NIMH) Data Archive (NDA, Collection #2846). This study was supported by a National Science Foundation CAREER Award (to A.P.M. under Grant Number 2045095), National Institute on Drug Abuse Grant 1R34DA050297-01 (A.P.M.), Jacobs Foundation Early Career Award (A.P.M.), CIFAR Global Azrieli Fellowship (to A.P.M.), Quad Fellowship (to M.N.), Nakajima Foundation Scholarship (to M.N.), National Science Foundation Graduate Research Fellowships (to A.L.B., U.A.T., and C.L.M. under Grant Number DGE-1845298), Behavioral and Cognitive Neuroscience Training Grant (to A.T.P. under Grant Number NIH T32-MH017168), as well as startup funds from the University of Pennsylvania. This manuscript reflects the views of the authors and may not reflect the opinions or views of the NIMH or of the submitters contributing original data to NDA.
The authors declare no competing financial interests.
This paper contains supplemental material available at: https://doi.org/10.1523/JNEUROSCI.0025-25.2025
- Correspondence should be addressed to Allyson P. Mackey at mackeya{at}upenn.edu.
