Abstract
Early structural and molecular development of the human cortex is extensively studied, but little is known about the development of neuronal activity across cortical regions. We used dense array electroencephalography recordings and a machine learning-based measure, functional brain age (FBA), to study spatiotemporally resolved maturation of cortical activity across the birth transition in human infants (male and female). We found clear spatial FBA gradients indicating more mature frontal cortical activity relative to other brain regions (geometric axis), as well as more mature activity in association cortices relative to sensory cortices (hierarchical axis). The frontal advance was explained by more mature bursting characteristics, a hallmark of early endogenous neuronal activity. The findings jointly support an advanced maturation of neuronal ensemble activity in cortical regions that are preparing to host synergistic, large-scale network interactions, a key global characteristic of mature brain function.
Significance Statement
Orderly brain organization from late pregnancy to early infancy is of paramount importance for optimal lifelong health. These complex processes are driven by brain activity, but surprisingly little is known about the maturation of this activity across brain regions. We used high-density electroencephalography recordings and machine learning-based measures of brain age to study regionally specific maturation in the developing brain. We showed that frontal and central regions are more mature than other brain regions with these differences largely driven by characteristics of bursts of brain activity. While challenging common ideas about delayed frontal development, our findings suggest that the early maturation in bursts of activity in frontal regions reflects its upcoming roles in the large-scale organization of the brain.
Introduction
Orderly neurodevelopment across the birth transition, from the last months of pregnancy to the first month after term age, is of paramount importance for lifelong neurocognitive success. Early brain organization relies on genetically programmed spatiotemporal passages of molecular expressions and cellular growth (Molnár et al., 2020; Bhaduri et al., 2021), which require concurrent neuronal activity for organizing vastly complex brain structures and functions (Martini et al., 2021; Wu et al., 2024). Therefore, characterizing arealization in the maturation of neuronal activity is fundamental to understanding the mechanisms that drive and support early activity-dependent brain organization, as well as regionally specific sensitivities to environmental or medical factors (Tokariev et al., 2019a; Taoudi-Benchekroun et al., 2022; Yrjölä et al., 2022; França et al., 2024).
Despite large strides in other directions of developmental neuroscience, very little is known about the early arealization in cortical neuronal activity (Molnár et al., 2020; Luhmann et al., 2022). The experimental literature focuses systematically on selected brain regions (Ackman et al., 2012; Mizuno et al., 2018; Suárez et al., 2023), precluding whole-brain regional comparisons. In turn, clinical studies using scalp-recorded electroencephalography (EEG) in prematurely born human infants have established weekly progression of cortical neuronal activity from the final trimester of pregnancy to the early weeks following term age (Parmelee et al., 1968; André et al., 2010; Wallois et al., 2021). The extensive neonatal EEG literature describes regional differences in the development of human cortical activity that align in some respects with arealizations reported in molecular (Bhaduri et al., 2021; Vogel et al., 2024), microscopic (Conel, 1939; Huttenlocher, 1990), and macroscopic (Hüppi et al., 1998; Batalle et al., 2017; Ball et al., 2020; Larivière et al., 2020; Vo Van et al., 2022) studies. However, studies on spatial ordering in the maturation of neuronal activity have been previously compromised by the lack of commensurate measures for all cortical regions, such as the recently introduced machine learning-based models for estimating functional brain age (FBA) from scalp-recorded EEG signals (Stevenson et al., 2020; Ansari et al., 2024).
Here, we estimated FBA at high spatial resolution to study arealization in the early maturation of cortical activity in human infants. We hypothesized that measures of FBA show clear spatial patterning with spatial gradients along geometric (Tokariev et al., 2019a; Pang et al., 2023) and/or hierarchical axes (Sydnor et al., 2021; Larsen et al., 2023), which frame global brain function throughout later life. Since neonatal EEG maturation is based on the time of conception rather than the time of birth, we studied EEG maturation in preterm-born infants over a period of 3 months across the “birth transition” (Ji et al., 2024), from 2 months before [∼32 weeks postmenstrual age (PMA)] to 1 month after a normal term birth (∼44 weeks PMA). By leveraging sequential high-density EEG recordings, source reconstruction methods, and computational FBA estimates, our study identified distinct spatial gradients in the maturation of neuronal activity across cortical regions that persist across the birth transition.
Materials and Methods
An overview of the study design is shown in Figure 1. We used high-density EEG recordings from 136 preterm infants measured at two time points between 33 and 44 weeks PMA. We converted EEG recordings (sensor space) into a source space representing 44 parcels of the cortical surface. A machine learning-based prediction of age (FBA) was estimated from each parcel (training) and then applied to the remaining parcels (testing). Differences between FBA on training and testing parcels were used to define relative maturation. This maturation metric was then reordered according to four axes (three geometrical and one hierarchical), and correlation analysis was performed to determine the presence or absence of linear gradients.
An overview of the study pipeline for assessing spatial gradients in maturation of cortical activity in human infants. A, Preterm-born infants were recorded at two timepoints using 124 channels of EEG, and the data were transformed to cortical source activities in 44 parcels. B, Algorithms for functional brain age (FBA) were trained for the parcel signals using 20 quantitative EEG features and support vector regression with PMA as the target. Relative maturity in cortical activity was measured as FBA difference (ΔFBA) between cortical parcels by training the FBA on one parcel, followed by testing it on all other parcels. C, This yielded a matrix of ΔFBA that depicts the pairwise FBA differences between each training and testing parcel; this representational image shows ΔFBA sorted along a frontal–occipital axis, and note that the diagonal elements are zero. D, The relative maturity was ordered according to the axis of interest, and correlation analysis was used to assess the presence of spatial gradients. E, Several biological and technical factors were examined as potential confounders of spatial gradients.
Description of the EEG dataset
The dataset was collected as part of a prospective, randomized controlled nursing trial (Welch et al., 2014) where previously published analyses indicated no global differences in FBA estimates between trial arms (Yrjölä et al., 2022); hence, all trial participants were pooled together in the present study. We received permission from the primary investigators of this clinical trial to use the EEG recordings in our study; the dataset is not publicly available. This randomized controlled trial was registered at ClinicalTrials.gov (#NCT01439269), and the study site at Columbia University Medical Center had approved all recruitment, consent, and study procedures (Welch et al., 2014; Yrjölä et al., 2022). Written informed consent was obtained from mothers before the start of intervention procedures. As previously described, each preterm infant (born at between 26 and 34 weeks gestational age) was recorded at least twice between 33 and 44 weeks PMA (37–42 weeks is considered term-equivalent age; Welch et al., 2014). The EEG was recorded from sleeping infants at least 30 min after feeding and between 11 A.M. and 4 P.M. using an Electrical Geodesics 200, 16-bit portable recording system, and a 128-electrode net (electrode impedances were maintained at the level recommended by the manufacturer during recording; Electrical Geodesics). The average recording duration was 62 min ± 2.4 min. The reference electrode was placed on the vertex, and signals were sampled at 1 kHz with bandpass filtering at 0.1–400 Hz. Analysis was performed, offline, on an average montage after excluding electrodes located on the face (124 channels). Here, our primary interest was an infant's PMA (in weeks) at the time of EEG recordings, while other demographics from the original trial were assessed as potential confounders of the FBA estimate: trial arm, biological sex (male or female), and neurodevelopmental outcome at 18 months of age (Bayley Scales of Infant and Toddler Development III).
EEG preprocessing
After exporting to European Data Format (EDF), all further processing was performed using custom scripts (Matlab Version 2020b, MathWorks). The EEG signals were initially filtered (Butterworth, fourth-order bandpass filter with a lower cutoff frequency of 0.5 Hz and an upper cutoff frequency of 32 Hz) and resampled to 256 Hz ensuring coverage of all frequencies seen in scalp-recorded preterm EEG activity. A process of automated artifact detection was applied with outliers of EEG amplitude and frequency content detected. EEG recordings with an excessive burden of artifact were excluded from further analysis. ICA-based artifact removal was then applied to the remaining recordings to remove noncortical components (Winkler et al., 2011): We automatically selected independent components for removal that were highly oscillatory and, therefore, associated with heart rate or respiration based on similar amplitude and frequency domain criteria as used in the artifact detection stage. EEG recording rejection, artifact detection, and removal were all performed before the training and testing of the FBA. The full details of artifact detection and removal are given in Text S1 and on GitHub (https://github.com/nstevensonUH/spatial_maturation_preterm_brain).
Reconstruction of cortical signals
EEG electrodes measure the electrical activity of a mixture of cortical sources across a wide spatial extent (Schoffelen and Gross, 2009; Schaworonkow and Nikulin, 2022). To improve the spatial specificity of our analysis, we used source reconstruction methods applied to the preprocessed scalp EEG data to estimate cortical parcel signals. We used the method described in Tokariev et al. (2019b) which includes a realistic three-shell (scalp, outer/inner skull) infant head model and dynamic statistical parametric mapping (Dale et al., 2000). The source space was represented by 8,014 dipoles uniformly distributed across the cortical surface and oriented perpendicular to it. All sources were then clustered into 58 parcels of near equal size (the number of sources per parcel was 129 ± 25, mean ± SD) that were symmetric across hemispheres. Parcel signals were computed as weighted means of source signals belonging to it (Tokariev et al., 2019b). The initial parcellation was reduced to 44 symmetric parcels by removing fourteen parcels with compromised fidelity of source-level signal reconstruction due to their distance from the scalp (Fig. S1).
Feature selection for FBA estimation
Computational EEG summary features were calculated on the entire EEG recording. We used the full recording duration as FBA estimation has been shown to be robust to, if not improved by, a mixture of sleep/vigilance states (Stevenson et al., 2020). The initial feature set included 49 features selected from the set used for training a global FBA estimate of preterm EEG recordings (Stevenson et al., 2020). The core of the feature set consisted of summary statistics of bursts that were identified by applying a data-driven, adaptive threshold to the EEG signal envelope (the absolute value of the analytic associate of the EEG (Roberts et al., 2014; see Fig. S2 for a visual guide on their calculation). These measures were supplemented with other measures of EEG activity from the analysis of amplitude, frequency, and information/complexity domains. Furthermore, features that did not contribute to age prediction (based on cross-validation analysis of predictions from each parcel) were removed, resulting in a final set of 20 features (Table 1).
The subset of EEG features used to predict age per parcel signal
Calculation of FBA
A combination of features generated an FBA for each cortical parcel (Stevenson et al., 2020). Training of support vector regression (SVR) models included data augmentation (synthetic data) to improve homogeneity in PMA distributions. Synthetic data were drawn from a multivariate normal distribution with mean and covariance estimated from the existing dataset with adjustments for age-dependent changes in features (adjustments via polynomial fitting with order selection via Akaike Information Criterion). For the SVR, the Matlab function fitrsvm was used with a coarse Gaussian kernel (kernel scale = 16, box constraint = interquartile range of PMA divided by 1.349, epsilon = interquartile range of PMA divided by 13.49, features were normalized). Further details on the training process are available in Text S1.
Analysis of spatiotemporal patterning in cortical functional maturation
Parcel-level heterochronicity
A matrix of FBA estimates was generated by training on all available data from one cortical parcel and testing on data from another parcel (when the training and testing parcel was the same, fivefold cross-validation was used to estimate the FBA). Differences in cortical maturation were assessed by comparing training FBAs to testing FBAs per parcel per recording (i.e., paired data), yielding a fully counterbalanced, unbiased interparcel comparison scheme.
Global-level FBA patterning
We studied two alternative hypotheses in global patterning: First, “geometric patterning” was motivated by recent studies in adults (Pang et al., 2023) and infants (Tokariev et al., 2019a) showing that brain geometry explains functional constellations in neuronal activity. This was studied with a hypothesis that the observed FBA gradients link to three simple, geometrical gradients: frontal–occipital (FO), lateral, and superior–inferior, respectively. Second, “hierarchical patterning” was motivated by recent studies showing that developmental gradients across the sensorimotor-to-association (S–A) axis are seen at multiple levels, from genetic expression to in vivo neuroimaging (Sydnor et al., 2021). Therefore, the S–A axis provides a physiologically relevant framework to study the cascades of functional neurodevelopment, including maturation of cortical neuronal activity. The S–A axis was derived by first annotating each parcel as belonging to either the sensorimotor (S) or association (A) region according to Larsen et al. (2023) and Ball et al. (2020). Each S or A parcel was ranked by aligning the S–A rank in Sydnor et al. (2021) with our parcellation (Fig. S3).
Computation of FBA gradients and statistical analysis
The underlying measure of maturational gradients was the ΔFBA per parcel. For a given parcel, an age prediction (FBA) was trained and then applied to all remaining testing parcels. ΔFBA was defined as the difference between all testing parcel FBAs and the training parcel FBA (evaluated in a cross-validation) averaged across all recordings and testing parcels. This resulted in an ΔFBA per parcel relative to all other parcels. To determine spatial gradients, ΔFBA was ordered according to the axis of interest, and a correlation analysis was performed to determine if there was a linear relationship between ΔFBA and axial distance (i.e., distance between parcels along a particular axis). Pearson's linear correlation coefficient was used with 95% confidence intervals estimated using a bootstrap resampling procedure. We also observe the line of best fit between ΔFBA and axial distance on a per recording basis; a predominance of positive or negative slopes in the linear fit supports groupwise findings. Differences in the correlation coefficients (Δr) between ΔFBA and axes were assessed using Zou's estimates of confidence intervals (level of significance was 0.05, overlapping correlations; Zou, 2007). A confidence interval of Δr that did not include 0 was deemed significant.
The use of FBA was supported with an assessment of FBA accuracy on both training (via 10-fold cross-validation) and testing parcels. Accuracy was defined using the mean absolute error between FBA and PMA and the correlation between FBA and PMA.
Analysis of potential confounders
The effects on spatial gradients were tested for both technical and biological/clinical confounders. Technical confounders relate to the process of estimating the FBA and spatial differences across the brain (notably, source space reconstruction), while biological/clinical confounders relate to biological differences within the cohort. Biological confounders were assessed by ensuring that a comparable, significant spatial gradient was found in each subgroup (biological sex, trial arm, neurodevelopmental outcome). Technical confounders were assessed using a null model of random spatial mixing of EEG channels to determine if the observed gradients could be explained by source space reconstruction. The null model contains the temporal variation of the EEG with PMA while mixing spatial structure in the sensor space.
Results
Accuracies in local FBA assessments are globally comparable
Age estimates were calculated for 44 cortical parcels in 267 recordings from 135 infants between 33 and 44 weeks PMA. The accuracy of age prediction within a parcel was high with a median MAE of 1.18 weeks (IQR: 1.13–1.22; n = 44; cross-validation results) and a median correlation between predicted age and PMA of 0.813 (IQR: 0.797–0.833; n = 44). The accuracy of FBA estimates was comparable across cortical areas (median MAE = 1.40 weeks; IQR: 1.26–1.64 weeks and median r = 0.784; IQR: 0.750–0.808, n = 1,892); thus, differences in FBA between cortical areas were assumed to represent genuine maturational offsets in neuronal activity.
Local FBA estimates show significant brain-wide gradients
Next, we studied whether spatial heterogeneity in the maturation of neuronal activity was ordered according to geometrical or hierarchical axes. Geometrically, we found a significant maturational gradient along the FO axis (r = 0.606, 95% CI: 0.377–0.765; n = 44; Fig. 2A) and lateral axis (r = −0.487, 95% CI: −0.685 to −0.226; n = 44; Fig. 2B). The FBA estimated from the frontal parcels was significantly more mature than the occipital FBA with a difference of 1.39 weeks (95% CI: 0.96–1.83 weeks) between frontal and occipital extrema. Along the lateral axis, central neuronal activity was more mature than temporal activity with a difference of 0.96 weeks (95% CI: 0.43–1.50 weeks; n = 44) between central and temporal extrema. However, no significant gradients were observed along the superior–inferior axis (correlation analysis; r = 0.246, 95% CI: −0.055 to 0.506; n = 44). These gradients were evidenced in individual recordings with consistent slopes seen in both the FO and lateral axes (Fig. 3). Hierarchically, there was a significant gradient of maturation along the sensorimotor–association (SA) axis (r = 0.398; 95% CI: 0.115–0.622; n = 44; Fig. 4A) and significant groupwise differences in maturation between parcels assigned as association (A) versus sensorimotor (S) cortex (median association cortex ΔFBA = 0.07 weeks; IQR: −0.47 to 0.41 vs median sensorimotor cortex ΔFBA = −0.50 weeks; IQR: −0.87 to −0.07; n = 44, p = 0.046; Fig. 4B). There was a moderate effect size in ΔFBA between sensorimotor and association parcels (Cohen's D = 0.77; 95% CI: 0.51–1.02). The correlation between ΔFBA and the FO axis was significantly higher than the correlation between ΔFBA and the SA axis (Δr 95% CI: 0.017–0.425, n = 44) even though FO distance and SA rank were correlated (r = 0.717; 95% CI: 0.534–0.836, n = 44). There was no significant difference between ΔFBA correlations with FO versus lateral axis (Δr 95% CI: −0.208 to 0.449, n = 44).
Gradients of cortical activity maturation with respect to geometric patterning. Relative maturation along (A) the frontal–occipital axis—zero is the vertex plane near the central sulcus, positive is frontal, and negative is occipital—and (B) the lateral axis; zero represents the sagittal midline; the higher the value, the more temporal the location (correlation is estimated for the absolute deviance from the midline to allow assessment of symmetry). C, The relative activity maturation (ΔFBA) in all parcels projected onto a 3D cortical surface model. In (B) and (C), the detrended data for all testing regions are shown in the background. The two filled circles in (A) define ΔFBA at exemplar regions used in Figure 6. In all figures, positive ΔFBA means more mature local cortical activity, and negative ΔFBA means less mature local cortical activity compared with all other parcels on average (F, frontal; L, left). Note that |ΔFBA| greater than 0.17 weeks corresponds to a p < 0.05 under a null hypothesis of |ΔFBA| = 0 based on the Wilcoxon signed rank test and corrected for multiple comparisons (Bonferroni's method).
Spatial gradients per recording. A, Frontal–occipital axis (94% of recordings have a positive slope), (B) lateral axis (93% of recordings have a negative slope), and (C) inferior–superior axis (75% of recordings have a positive slope). A total of 231 recordings had outlier-free data in more than 22 parcels (out of 44) and were included in this analysis. Slope is measured in weeks per centimeter.
Gradients of cortical activity maturation with respect to hierarchical patterning. A, Relative maturation along the sensorimotor–association (SA) axis where S on the x-axis denotes the highest rank of the sensorimotor axis, A on the x-axis denotes the highest rank of the association axis, and 0 is the transition at which parcels were assigned as either sensorimotor or association. B, Groupwise differences in ΔFBA between sensorimotor (S) and association (A) parcels (asterisk denotes p < 0.05; Mann–Whitney U test). C, The spatial distribution of SA ranking used in (A) and (B) projected on a 3D cortical surface model (FL, front left).
Spatial FBA gradients persist throughout the birth transition
Maturational gradients were consistent throughout the studied period of 2–3 months across the birth transition. The correlations between ΔFBA and axial distances were comparable in infants measured before 36 weeks PMA and after 36 weeks PMA, respectively (Table 2 and Fig. S4). Potential biological confounders did not affect the spatial distribution of FBA and, therefore, maturational gradients along any of the four axes tested (Table S1). Likewise, spatial FBA gradients were not affected by the technical confounders related to the process of source space reconstruction (Table S2).
Gradients estimated using ΔFBA calculated on all available data, only recordings from before 36 weeks PMA, and only recordings after 36 weeks PMA
Maturation of neuronal bursting drives the FBA gradients
Examining the computational EEG features used to train FBA models reveals potential mechanistic underpinnings of the observed FBA gradients. To explain the FBA gradients, the given feature should correlate significantly with both the infant's age and the given axial distance, unaffected by corrections for multiple comparisons (Fig. 5). Only two features were correlated significantly with both the FO axis and an infant's age, independently: sharpness of the bursts with 1–2 s duration (Feature 10) and slope of the line of best fit applied to burst duration (BD) versus burst area (Feature 11). In addition, two features were associated with lateral distance and age, independently: sharpness of bursts with 125–250 ms duration (Feature 9) and intercept of the line of best fit applied to burst duration versus kurtosis (Feature 14). Only a single feature was associated with SA rank and age, independently (Feature 10). The use of machine learning-based combinatorial measures more accurately represents age in neuronal activity and, therefore, maturation providing a more accurate evaluation of spatial patterning (Fig. 6).
Linking neuronal activity to ΔFBA gradients. A, Correlations between EEG features (x-axis) and distance along the FO axis, ordered according to |r|. B, Correlations between EEG features (x-axis) and infant age. To provide a mechanistic explanation for the observed spatial FBA gradients (ΔFBA), a feature must show statistically significant correlation to both the given spatial patterning (A) and the infant's age (B). The marker colours denote feature classes. The exact descriptions of the feature IDs are detailed in Table S2. The horizontal dashed line denotes the threshold for statistical significance (p = 0.05; corrected for multiple comparisons in the case of individual feature analysis).
Improving the evaluation of spatial differences in maturation using a combinatorial measurement of age (FBA) compared with underlying individual features. A, Top, FBA trained on an exemplar frontal parcel and tested on an exemplar occipital parcel (F) versus FBA trained on an exemplar occipital parcel and tested on an exemplar frontal parcel (O). Bottom, the distribution of FBA around the line of best fit with the frontal parcel showing more mature activity (Cohen's D between parcels is 3.03, 95% CI: 2.75–3.31). B, Burst sharpness (BS; in bursts of 1–2 s in duration) in an exemplar frontal parcel (F) and an exemplar occipital parcel (O). Bottom, The distribution of BS around the line of best fit with the frontal parcel showing more mature activity (Cohen's D between parcels is −2.03, 95% CI: −2.23 to −1.82). C, The slope of BS versus burst duration (BD) in an exemplar frontal parcel (F) and an exemplar occipital parcel (O). Bottom, the distribution of the slope around the line of best fit with the frontal parcel showing more mature activity (Cohen's D between parcels is 0.36, 95% CI: 0.19–0.53). The location of the exemplar parcels is shown in Figure 2.
Discussion
Our findings indicate clear spatial gradients in the early development of human cortical activity, with frontal activity more mature than occipital regions by approximately 1 week. Prior literature has reported spatiotemporal patterning in early human brain development ranging from gene expression, micro- and macrostructural measures, and global connectomics (Huttenlocher, 1990; Hüppi et al., 1998; Bhaduri et al., 2021; Vo Van et al., 2022). Our present work extends prior literature by describing brain-wide interareal differences in the maturation of neuronal activity, the driver of activity-dependent brain organization during late pregnancy and early infancy (Molnár et al., 2020). A key strength of our study is the use of source space reconstruction from high-density EEG with more electrodes than cortical parcels (overdetermined spatial sampling) to ensure spatial fidelity was maximized and machine learning to ensure characteristics most associated with age were extracted from cortical signals.
The present findings of advanced frontal development conflict with the long history of neuroscience literature showing protracted postnatal development of the human frontal lobe structures and functions (Thatcher, 1991; Kostović et al., 2019; Kolk and Rakic, 2022). Our present work differs from the prior literature in two key ways, we studied: (1) the activity of the cortical neuron ensembles, rather than structural organization or molecular expressions, and (2) a period of birth transition, rather than a later postnatal period which is characterized by the late emergence the higher cognitive functions. It is, therefore, entirely possible that the neuronal ensembles in the frontal cortex exhibit activity that is mature with respect to other regions of the brain while showing considerably protracted development of higher cognitive qualities. Our findings reinforce the concept of high spatiotemporal heterogeneity in the developing brain, with apparently contrasting trajectories of growth in microstructure, macroscopic measures, and different molecular markers (Hüppi et al., 1998; Ball et al., 2013; Bhaduri et al., 2021; Sydnor et al., 2021).
Global-level FBA gradients find putative parallels in two recently introduced frameworks: geometric and hierarchical. Geometrically, the frontal–occipital and lateral gradients resemble anatomy-derived spatial eigenmodes which shape large-scale functional connectivity and spatial brain modes in both adults (Pang et al., 2023) and infants (Tokariev et al., 2019a). Frontal neuronal activity in the human newborn hosts the first synergistic scaffolds (Varley et al., 2025) and shows distinct brain-wide connectivity with developmental sensitivity to medical or environmental effects (Videman et al., 2017; Taoudi-Benchekroun et al., 2022; Yrjölä et al., 2022). Moreover, a large number of structural (Kersbergen et al., 2014; Batalle et al., 2017) and molecular (Bhaduri et al., 2021; Vogel et al., 2024) studies have repeatedly indicated the existence of measure-specific developmental axes where frontal regions stand out either alone (Ji et al., 2024) or together with temporal and parietal regions (Taoudi-Benchekroun et al., 2022). Hierarchically, the observed FBA gradients largely conform to a systems-level organization where the principal axis divides the cortex into sensory and higher-order association regions (Ball et al., 2020; Sydnor et al., 2021; Larsen et al., 2023). Recent studies have established evidence for distinct developmental trajectories in the sensory and association systems, respectively, with neurodevelopmental cascades of critical periods culminating in the cortical association areas (Cadwell et al., 2019; Larsen et al., 2023).
Strict statistical comparisons favor a geometric over a hierarchical axis although both are clearly compatible with the present data. The relatively weaker relationship with the SA axis appears to result from relatively lower ΔFBA in the parietal and inferior temporal areas of the association cortex; however, parcel-wise FBA estimates in these cortical regions may represent mixed effects due to the rich functional diversity in these subregions (Krubitzer and Prescott, 2018; Petersen et al., 2024). From a functional perspective, the hierarchical SA axis would offer a more powerful framework to reconcile our FBA findings with a larger body of neurophysiological literature: Early EEG-measured cortical activity emerges first in the sensory cortices where peripheral sensory inputs trigger local bursting (neuronal ensemble activity) to train local cortical information processing (Cadwell et al., 2019; Molnár et al., 2020; Luhmann et al., 2022). In contrast, frontal cortical activity measured with EEG emerges relatively late (Wallois et al., 2021) with inputs to frontal regions primarily from other cortices and subcortical structures with only weak modulation from sensory systems (Gómez et al., 2023). Recent cortex-wide calcium imaging studies have also indicated a clear global topography that complies with the SA axis (Warm et al., 2025). Therefore, the cortical activity of the frontal regions develops directly to host synergistic ensembles that support global brain function throughout life (Luppi et al., 2024; Varley et al., 2025) while altered frontal neuronal activity reflects mental and cognitive compromise (Chini and Hanganu-Opatz, 2021). Nevertheless, the geometric interpretation of maturation is equally valid as geometric constraints (e.g., exponentially decaying wiring rules) have been shown to parsimoniously explain numerous features of large-scale structural (Ercsey-Ravasz et al., 2013; Henderson and Robinson, 2013; Markov et al., 2013; Roberts et al., 2016) and functional connectivity (Robinson et al., 2016; Pang et al., 2023). In particular, cortical activity measured with EEG is dominated by long-wavelength modes (Nunez, 2000; Tokariev et al., 2019a) which increase the robustness of neural field descriptions as they are relatively insensitive to finer-scale cortical variations (Gabay and Robinson, 2017).
Regarding the underpinning neural mechanisms, the objective machine learning-based FBA models were found to favor nontraditional EEG features selecting alternative characteristics of EEG bursts, the hallmark of the early endogenous cortical activity (Vanhatalo and Kaila, 2006; Molnár et al., 2020; Wallois et al., 2021), which shows continuous and monotonic changes with age (Stevenson et al., 2020). The strongest EEG feature was the measure of burst sharpness in both the geometrical and hierarchical axes. Studies modeling cortical activity in silico as a stochastic process have implied that increased burst sharpness relates to decreased damping within the burst-generating mechanism and, therefore, increased activity/excitability in the cortex (Roberts et al., 2014). A more physiologically detailed model has also implicated metabolic and ionic resources as modulators of burst sharpness (Dutta et al., 2023).
Such burst characteristics offer an optimal temporal frame for brain-wide integrative processing in frontal areas. Comparable characteristics would emerge later in the sensory areas after they transition from robust local sensory responses (Molnár et al., 2020) to processing wider corticocortical integration (Petersen et al., 2024). Drawing from this rationale, our findings support a mechanistic theory where upcoming functional integrative roles of the frontal and other association regions are reflected in the more mature activity of local neuronal ensembles (higher FBA). Such developmental cascades are also seen in frontal cortical activity networks, where developmental advance or compromise links selectively to upcoming neurocognitive performance (Tokariev et al., 2019a; Bitzenhofer et al., 2021; Chini and Hanganu-Opatz, 2021; Luppi et al., 2024; Yrjölä et al., 2024). Conversely, the relatively lagging maturation in the sensory cortices would reflect the longer gradual shifting from the fetal mode of sensory-driven activity to cortical endogenous activity (Molnár et al., 2020; Wallois et al., 2021). Although this hypothesis is essentially circumstantial, it provides a testable framework for future experimental and modeling studies.
Our findings of spatial gradients in functional brain development across the birth transition must be tempered by the fact that premature birth is not normal, and early exposure to the ex utero environment subtly affects the EEG (Scher et al., 1994). We assume that these changes do not alter the existence and direction of the spatial gradients of maturation found in our study. We also chose to include all infants in the cohort to maximize the data size for training machine learning algorithms; ideally, only infants with typical neurodevelopment through this period should be included. We showed that potential confounders, such as sex, outcome, and intervention, did not affect our findings, most likely due to the machine learning methods adjusting for confounders during training in pursuit of maximum age prediction accuracy (EEG characteristics are associated with confounders; Griesmaier et al., 2014; Stevenson et al., 2020). While our study focused on detecting spatial gradients of maturation in the developing brain, potential disruption to these gradients from injury or disease offers an intriguing avenue of investigation for a prospective study. Any additional data collection also provides the opportunity to evaluate the reproducibility of our findings across alternative populations.
Footnotes
We thank the team of Profs. Marta Welch and Michael Myers at Columbia University, USA, for providing access to their EEG recordings from the Family Nurture Intervention trial. This work was supported by the National Health and Medical Research Council of Australia (#2002135), Academy of Finland, Sigrid Juselius Foundation, and Finnish Pediatric Foundation (Lastentautien Säätiö). The APC was funded by the Pediatric Research Center, HUH. The funders did not have any role in study design, data collection, data analyses and interpretation, or writing of the report.
The authors declare no competing financial interests.
This paper contains supplemental material available at: https://doi.org/10.1523/JNEUROSCI.0610-25.2025
- Correspondence should be addressed to Nathan J. Stevenson at nathan.stevenson{at}qimrb.edu.au.
