Regional BOLD variability reflects microstructural maturation and neuronal ensheathment in the preterm infant cortex

Domain: neuroscience
Published: 2026-04-09
Editorial outcome: major_revision

Abstract
Blood Oxygen Level Dependent (BOLD) variability reflects meaningful brain activity, yet its structural and biological correlates during early development remain unknown.
Using longitudinal resting-state fMRI and multi-shell diffusion imaging acquired longitudinally in 54 very preterm infants (at 33-weeks' gestational age and term-equivalent-age) and 24 full-term newborns, we investigated how BOLD variability evolves in very preterm infants, its relationship with cortical microstructure and gene expression, using the BrainSpan dataset, and how it differs from full-term newborns at term-equivalent age.
During preterm development, BOLD variability increased in primary sensory-sensorimotor and proto-Default-Mode-Network regions, accompanied by decreases in cortical diffusivity. Gene expression analysis revealed concurrent upregulation of genes mediating gliogenesis and neuronal ensheathment. At term-equivalent age, very preterm infants showed decreased BOLD variability and increased cortical diffusivity, compared to full-term newborns.
In this work, we show that BOLD variability reflects cortical microstructural maturation, mediated by upregulation of gliogenesis and neuronal ensheathment. Interruption of these processes by preterm birth identifies putative mechanisms of preterm brain injury.

Introduction
Neural dynamics can be measured in a non-invasive, though indirect way, using functional MRI (fMRI) through the Blood Oxygen Level Dependent (BOLD) signal, measured with functional MRI. Resting-state fMRI has been a particularly attractive paradigm to obtain richly structured spontaneous activity. So-called “static” analysis techniques obtain measures that summarize the average behavior over a whole run.
However, the “BOLD signal variability,” characterized by its variance, captures the magnitude of moment-to-moment regional variations in BOLD signal (Garrett et al., 2010, 2011, 2013). It is thought to indirectly reflect spontaneous brain activity arising from moment-to-moment fluctuations in neuronal activity. In fact, BOLD variability has been shown to correlate with electrophysiological measures of neural dynamics, including EEG power and temporal signal variability (Baracchini et al., 2021; McIntosh et al., 2010; Misic et al., 2011; Uddin, 2020; Vakorin et al., 2011). These temporal fluctuations in neuronal activity, as measured by electrophysiological and EEG studies, are thought to contribute to synaptic connectivity and associate with better cognitive performance (Faisal et al., 2008; McIntosh et al., 2008; Misic et al., 2011; Uddin, 2020; Vakorin et al., 2011), thereby underlying brain function.
Although BOLD variability is affected by vascular factors such as baseline cerebral blood flow and neurovascular coupling (Whittaker et al., 2016), studies combining BOLD and arterial spin labeling (ASL) imaging show that group differences remain even after accounting for these vascular influences (Garrett et al., 2017; Roberts R.P., 2023). This evidence indicates that BOLD variability reflects not only vascular effects but also neural dynamics, supporting its use as a meaningful index of underlying functional brain activity.
In adults, BOLD variability varies across cortical regions, with greater variability in association and transmodal regions (Baracchini et al., 2021). Greater BOLD variability has been linked to increased network organization, greater cognitive performance and increased task difficulty across different studies (Baracchini et al., 2021; Baracchini et al., 2023; Garrett et al., 2011, 2013; Goodman et al., 2024; Nomi et al., 2017; Protzner et al., 2013; Roberts R.P., 2023; Zhang et al., 2022), highlighting its functional relevance.
BOLD variability typically increases during childhood (Wang et al., 2021) and serves as a strong predictor of age across the lifespan (Garrett et al., 2010), then declining with aging (Baracchini et al., 2023; Garrett et al., 2013; Grady & Garrett, 2014; Rieck et al., 2022). To date, no studies have examined changes in BOLD variability in infants prior to term-equivalent age (TEA).
Additionally, few studies have explored the structural and biological underpinnings of BOLD variability, although evidence suggests it correlates with contiguous white matter (WM) microstructural integrity, both in adults and children (Burzynska et al., 2015; Good et al., 2020; Wang et al., 2021). Recently, using the ex-vivo BigBrain histological dataset (Amunts et al., 2013), Baracchini et al. found that BOLD variability was increased in cortical areas with greater laminar differentiation and higher neuronal density, particularly in granular cortical regions where the distinct layer IV supports efficient sensory processing and a broader range of functional responses (Baracchini et al., 2023). However, to date, no studies have examined the relationship between BOLD variability and intra-cortical microstructure, nor its biological correlates, specifically the underlying gene expression patterns, in adults or infants.
Diffusion-weighted imaging (DWI) measures water diffusion in brain tissue at the mesoscopic scale (Basser, 1995), enabling the study of brain microstructure. However, investigating gray matter (GM) microstructure with diffusion MRI (dMRI) poses significant challenges. Unlike white matter (WM), GM does not display orientational coherence, is not as heavily myelinated, and is characterized by a complex intra- and extracellular environment, with multiple barriers and different cellular compartments.
Traditional signal representation models of diffusion, such as Diffusion Tensor Imaging (DTI), provide insights into GM organization by measuring the anisotropic diffusion of water molecules within the tissue, without intrinsic diffusivity assumptions, capturing variations in diffusion patterns (Le Bihan et al., 2001). However, DTI assumes a Gaussian distribution of diffusion, which does not hold true in many biological tissues. To quantify the non-Gaussianity of diffusion, Diffusion Kurtosis Imaging (DKI) has been proposed as an extension of DTI, accounting for the non-Gaussian signal decay and estimating diffusion kurtosis as a reflective marker for tissue heterogeneity (Jensen et al., 2005). Studies have shown that diffusion kurtosis is more sensitive to microstructural complexity within tissues than DTI (Cheung et al., 2009; Paydar et al., 2014). Both DTI and DKI lack specificity for multiple diffusion signals within a single voxel (e.g., crossing fibers) and complex tissue microstructural features. To overcome this, multicompartment models have been proposed, to provide a more complete characterization of tissue complexity. The Spherical Mean Technique (SMT) model estimates microstructural parameters from the direction-averaged diffusion signal, hence removing directionality dependence, and provides measures sensitive to cellular composition, including intra‐axonal and extra‐axonal compartments (Kaden et al., 2016). Unlike other multi-compartment models, such as NODDI (Neurite Orientation Dispersion and Density Imaging), SMT does not rely on strong assumptions regarding fixed diffusivities, making it more suitable for complex tissue environments, like the GM (Zhang et al., 2012). Taken together, the combination of different dMRI models might provide valuable insights into GM microstructure and its association with BOLD variability.
While imaging can reveal the relationship between BOLD variability and cortical microstructure maturation during early brain development, gene expression analysis can provide mechanistic insight into the processes and pathways underlying developmental functional and structural changes occurring in the cortex. We hypothesized that regional changes in BOLD variability, coupled with changes in the cortical microstructure, are associated with the expression of genes linked to neocortical organization.
In addition to the lack of evidence about how BOLD variability changes during early brain development, few studies have examined the impact of preterm birth on its developmental trajectory (Amunts et al., 2013; Wu et al., 2016). Preterm birth is known to occur during a critical period of activity-dependent plasticity and brain development, characterized by axonal and dendritic growth and arborization, synaptogenesis, neural organization and myelination (Knickmeyer et al., 2008; Kostovic & Jovanov-Milosevic, 2006). Preterm birth exposes the infants to a dramatic environmental change and leads to altered structural and functional brain development with adverse, long-term outcomes (Dubois et al., 2014; Huang et al., 2006; Huppi et al., 1998; Lordier et al., 2019; Nossin-Manor et al., 2013; Sa de Almeida et al., 2023; Sa de Almeida et al., 2020; Volpe, 2001, 2009).
In this study, we aimed to investigate the maturation of BOLD variability in very preterm (VPT) infants, longitudinally, from 33 to 40-weeks gestational age (GA), across different resting-state networks (RSNs). In addition, using a comprehensive dMRI analysis, combining metrics from DTI, DKI, and SMT models, we evaluated how observed changes in BOLD variability across the cortex align with regional measures of cortical microstructure. Furthermore, to better understand the biological correlates of the observed functional and microstructural changes, we examined spatiotemporal patterns of gene expression, in postmortem tissue samples over the same time period. Finally, we aimed to assess whether preterm birth impacts the expected developmental maturation of BOLD variability and cortical microstructure, in comparison to full-term (FT) birth.

Results
Patient demographics
In both resting-state fMRI (RS-fMRI) and dMRI samples, significant differences between the VPT and FT groups were observed, as expected, in the following perinatal variables: GA at birth, birth weight, birth height, birth head circumference, APGAR at 1 and 5 minutes and incidence of bronchopulmonary dysplasia (BPD). There were no differences between groups in sex, GA at 2nd MRI scan (TEA), neonatal asphyxia, intrauterine growth restriction, intraventricular hemorrhage grade I, and socio-economic parental status score (Largo et al., 1989), (Table 1 and Supplementary Table S1). There were no significant differences, for any perinatal clinical variable, between the RS-fMRI and dMRI samples, in the VPT infants’ group or the FT newborns’ group (Table 1 and Supplementary Table S2).

Longitudinal cortical BOLD variability and diffusion microstructural changes during early preterm brain development
In the VPT cohort, BOLD variability (estimated as BOLD SD) increased significantly (false discovery rate (FDR)-corrected; p < 0.05) from 33- to 40-weeks GA (wGA) in the primary sensory (Visual, Auditory), Sensorimotor and proto-Default-Mode-Network (pDMN) (composed of the Posterior Cingulate Cortex and Precuneus) (group 1). In contrast, BOLD variability did not change significantly in Paralimbic, Thalamus, Limbic or Prefrontal Cortex networks (group 2), (Figure 1 and Supplementary Table S3).
Complementary analysis using ALFF and fALFF to estimate BOLD variability yielded consistent longitudinal patterns across the same regions (see Supplementary Figure S1).
We next examined changes in diffusion microstructural metrics in VPT infants, from 33- to 40-wGA, across all cortical regions. Significant increases in BOLD variability were associated with an overall decrease in cortical diffusivity and kurtosis metrics (group 1, Figure 2).
In group 1, where BOLD variability increased significantly, SMT metrics globally decreased, namely intrinsic diffusivity (diff), extra-neurite mean diffusivity (extraMD), extra-neurite transverse diffusivity (extraTrans) and intra-neurite volume fraction (intra). Additionally, DTI showed significant decreases of both fractional anisotropy (FA) and mean diffusivity (MD) in all these regions. DKI revealed significant decreases in radial kurtosis (RK) and/or mean kurtosis (MK) across these group 1 regions.
In group 2, where BOLD variability did not change significantly, the microstructural changes were less coherent. More specifically, the thalamus displayed a significant increase in FA, accompanied by significant decreases in diff, extraMD, extraTrans and MD between 33- and 40-wGA. In contrast, the PFC displayed a marked decrease in both FA and intra-neurite volume fraction, accompanied by an increase in MD. Both limbic and paralimbic regions displayed an overall decrease in SMT diffusivities with decreased intra-neurite volume fraction and FA (Figure 2).
We conducted an additional analysis, clustering cortical regions based on the similarity of their longitudinal microstructural maturation, rather than BOLD variability, to compare the microstructural clusters with those derived from BOLD changes. We found that the microstructural patterns are consistent with the functional partition (Supplementary Figure S3).

Relationship between longitudinal cortical microstructural changes and BOLD variability delta changes
To identify the microstructural metrics that best predicted the observed changes in BOLD variability, we employed Lasso-based sparse multivariate regression models on the regional-average data. Out of the 9 microstructural metrics, 4 features (MD, RK, AK and FA) were retained, at the optimal pure LASSO penalty (λ₁se = 0.015, α = 1), as predictors of regional BOLD variability changes from 33- to 40-wGA. MD showed the strongest association (β = 0.419), followed by RK (β = 0.402), AK (β = 0.353), and FA (β = 0.332). The model achieved a high training fit (R² = 0.965); however, leave-one-region-out cross-validation indicated that these longitudinal changes explained 28% of the observed variance in BOLD variability (R²ₗₒₒcᵥ = 0.28, RMSE = 0.054). Using less penalized models - a Lasso Elastic Net (α = 0.5, λ_min = 0.0032) and a - additional metrics (intra, diff, extraMD, MK) were retained, increasing the explained variance to 61–64% (R²ₗₒₒcᵥ = 0.61–0.64, RMSE = 0.0295–0.032, respectively).

Genetic expression patterns during early brain development
Using a developmental transcriptomic dataset of bulk tissue mRNA data sampled from cortical tissue in 18 post-mortem prenatal human specimens aged 8-37 post conceptional weeks (pcw) (http://development.psychencode.org/), we compared genetic expression patterns between group 1 (regions where BOLD variability increased significantly from 33- to 40-wGA) and group 2 (regions with no significant longitudinal BOLD variability changes), from mid- to late-fetal period. We hypothesized that regions with changes in BOLD variability (group 1), which also presented coupled cortical microstructural maturational changes, would display an increase in expression of genes linked to neocortical organization.
From mid- to late-fetal period (16-22 to 35-37 pcw), 786 out of the total 5,287 genes displayed a significant group*time interaction, after correction for multiple comparisons (FDR-corrected; p < 0.05). From these genes, using pairwise comparisons to filter those with differential expression between group 1 and 2 from mid- to late-fetal, we identified 132 genes that increased significantly more in group 1 vs 2, whereas only 9 genes increased significantly more in group 2 vs group 1. Additionally, 64 genes were found to decrease significantly more in group 1 vs 2, and 45 genes decreased significantly more in group 2 vs 1.
Using geneset enrichment analysis (Wang et al., 2017), with a background set of 5,287 genes (Ball et al., 2020), we identified significant enrichments for several important neocortical organizational processes in genes that increased significantly more in group 1 than 2. These included: “neuronal ensheathment” (FDR = 0.000014, enrichment = 6.15, 14 genes), “gliogenesis” (FDR = 0.018, enrichment = 2.85, 14 genes, from which 7 genes are common to neuronal ensheathment), “extracellular structure organization/external encapsulating structure organization” (FDR = 0.018, enrichment = 3.58, 10 genes), “cytokine-mediated signaling pathway” (FDR = 0.008, enrichment = 3.6, 13 genes), “coagulation/wound healing” (FDR = 0.008/0.01, enrichment = 4.09/3.03, 10/14 genes, where the 10 genes from coagulation are common to both categories), “extrinsic apoptotic signaling pathway” (FDR = 0.01, enrichment = 3.8, 10 genes) and “muscle cell proliferation (FDR = 0.02, enrichment = 3.8, 9 genes) (Figure 3c). A list of all significant categories and their assigned genes can be found in Supplementary Table S4. No significant enrichment analysis results were found for the genes that increased significantly more in group 2 than 1, or for the genes that decreased significantly differently between group 1 and 2.
The combined list of enriched genes associated with neuronal ensheathment and gliogenesis categories comprised a total of 27 genes, including constituents of the myelin sheath (MBP, MAG, MAL, MOBP, PLLP, PLP1, PMP22, CLDN11), genes with a role in myelination regulation (NKX6-2, UGT8, GPR17), important for oligodendrocyte function and expression (MBNL2, APOD, VWA1), for tissue integrity (COL6A1), with a role in growth factor signaling (ERBB3), neuroprotection and/or glioprotection (TNFRSF1B, GPR37L1, IL6ST), neuronal and glial differentiation and development (GPR183, APCDD1, NDRG2, CSRP1, TSC22D4, CD9) and in fatty acid and sphingolipid biosynthesis (ELOVL1, PADI2). ANOVA followed by Tukey’s test identified a significant increase in expression of this group of genes in both groups (1 and 2) specifically from mid- to late-fetal (T1 to T2, group 1: z-ratio = -9.84, p-value = 5.88 × 10⁻¹⁴, group 2: z-ratio = -7.15, p-value = 1.27 × 10⁻¹¹). No significant changes were observed from early- to mid-fetal in both groups (T0 to T1, z-ratio = -0.59 to -0.33, p-value > 0.05). The increase from T1 to T2 was significantly higher in group 1, compared to group 2 (z-ratio = -9.77, p-value < 1.5 × 10⁻²²). At the late-fetal time-point (T2), these same genes were also expressed at significantly higher levels in group 1 compared to group 2 (z-ratio = 7.74, p-value = 1.80 × 10⁻¹³). In contrast, there were no differences in the expression of these genes between groups at T0 (z-ratio = -2.7, p-value = 0.073) and at T1 (z-ratio = -2.23, p-value = 0.22), (Figure 3d). The expression of each of these genes at each time-point and per group is shown in Figure 3e.

Effect of preterm birth on BOLD variability and cortical microstructure at TEA
Compared to FT newborns, VPT at TEA showed decreased BOLD variability in the thalamus, PCC, visual, auditory and limbic networks (Figure 4a). Complementary analysis using ALFF and fALFF to estimate BOLD variability yielded similar cross-sectional patterns across the same regions (see Supplementary Figure S4). Regarding cortical microstructural differences, in comparison to FT, VPT at TEA displayed overall increased diffusivities (extra-neurite mean diffusivity, extra-neurite transverse diffusivity, intrinsic diffusivity and Mean Diffusivity) across all RSNs (Figure 4b-c). VPT at TEA also presented increased kurtosis (RK and/or MK) across RSNs cortical regions, apart from the thalamus (Figure 4d). There were no significant differences between groups regarding AK across all cortical regions. Details of all regional differences per diffusion metric, between VPT and FT infants, are provided in Supplementary Table S5.

Methods
Subjects and study design
Research Ethics Committee approval was granted by the “Cantonal Research Ethics Committee of Geneva” and hosted at University Hospitals of Geneva (HUG). Written parental consent was obtained prior to infant’s participation to the study.
54 VPT infants (<32 wGA at birth) and 24 FT infants were recruited at the neonatal and maternity units of the University Hospitals of Geneva (HUG), Switzerland, from 2017 to 2020. Exclusion criteria for all babies included major brain lesions detected on the MRI, such as high-grade intraventricular hemorrhage, leukomalacia, as well as micro or macrocephaly or congenital syndromes.
For FT newborns, inclusion criteria included birth after 37 wGA, height, weight and head circumference above the 5th and below the 95th percentiles, APGAR score > 8 at 5 minutes and absence of resuscitation, infection or admission to the NICU (neonatal intensive care unit).
Six VPT infants were excluded from the study due to parental withdrawal or detection of a genetic anomaly, and two FT infants were excluded due to parental withdrawal.
Out study comprises both longitudinal and cross-sectional designs. The longitudinal design involves only our VPT infant cohort, who underwent MRI examinations at two time-points, firstly during the 33rd week of GA, and secondly at TEA. The cross-sectional design involves both the VPT at TEA and the FT newborns, who underwent a single MRI examination few days after being born (term age).
Infants whose MRI protocol acquisition was incomplete (not comprising a T2-weighted image, resting-state fMRI (RS-fMRI) sequence and multi-shell diffusion imaging (MSDI) sequence), without both longitudinal time-points (in preterm infants’ case) or whose images presented excessive motion were excluded from the analysis. The final sample comprised 46 VPT infants and 21 FT newborns, of which 31 VPT and 19 FT newborns had data suitable for the RS-fMRI analysis, while 39 VPT infants and 18 FT newborns had data suitable for the diffusion microstructural analysis, with an overlap of 24 VPT infants and 16 FT newborns between the two analyses. The flow chart of the participant selection process is provided in Supplementary Figure S5.

MRI acquisition
All infants were scanned after receiving breast or formula feeding during natural sleep (no sedation used). At 33 wGA, preterm infants were scanned using an MR compatible incubator (Lammers Medical Technology, Lübeck, Germany) and monitored using a Philips MR patient monitor Expression MR400 and Philips quadtrode MRI compatible neonatal ECG electrodes. At TEA, infants were scanned using a vacuum mattress for immobilization. All infants, at both time points, were monitored for heart rate and oxygen saturation during the entire scanning time. MR-compatible headphones (MR confon, Magdeburg, Germany) were used to protect the infants from scanner noise.
MRI acquisition at both time-points was performed on a 3.0T Siemens Magnetom MR scanner (Siemens, Erlangen, Germany), using a 16-channel neonatal head coil. T2-weighted images were acquired using the following parameters: 113 coronal slices, TR= 4990 ms, TE= 160 ms, flip angle= 150°, matrix size= 256 × 164; voxel size= 0.8 × 0.8 × 1.2 mm3. RS-fMRI data acquisition was obtained by means of T2*-weighted gradient echo echo-planar imaging (EPI) sequence with the following parameters: 590 images, TR = 700 ms, TE = 30 ms, 36 slices, voxel size = 2.5 x 2.5 x 2.5 mm3, flip angle = 60°, multiband factor = 4. Multi-shell diffusion imaging (MSDI) was acquired with a single-shot spin echo echo-planar imaging (SE-EPI) Stejksal-Tanner sequence with the following parameters TE=85 ms, TR=3170 ms, voxel size 1.8 ×1.8 ×1.8 mm3, multi-band factor = 2, GRAPPA 2. Images were acquired in the axial plane, in anterior-posterior (AP) phase encoding (PE) direction, with 4 volumes without diffusion-weighting (b0); 10 non-collinear directions with b= 200 s/mm2, 30 non-collinear directions with b= 1000 s/mm2; 50 non-collinear directions with b= 2000 s/mm2. Additional b0 images were collected with reversed phase-encode blips, AP and posterior-anterior (PA), resulting in pairs of images with distortions going in opposite directions.

RS-fMRI analysis pre-processing
The RS-fMRI data were preprocessed using SPM12 (Wellcome Department of Imaging Neuroscience, University College London, United Kingdom). The preprocessing steps included realignment to the mean functional volume, adjusting for motion, co-registration to the time-point specific structural image, alignment in Montreal Neurological Institute (MNI) space, and normalization, in which all the scans are warped to our cohort 40 wGA template. Slice-timing correction was omitted due to the fast TR of 700 ms used. Finally, the data were smoothed with a Gaussian filter of full width at half maximum (FWHM) of 5 mm. All volumes with a frame-wise displacement (Power, Mitra et al. 2014) greater than 0.5 mm or with a rate of BOLD signal changes across the entire brain (DVARS) greater than 3% were removed, along with the two previous and the two subsequent images. The remaining images were included for further analysis. A minimum of 50% of volumes remaining was set as sufficient criteria for inclusion.

Resting-state atlas creation
Resting-state networks (RSNs) were generated from our cohort of infants, comprising both VPT at 33- and 40 wGA, as well as the FT infants, and combined into a brain cortical parcellation atlas. To extract independent spatial networks, a group-level independent component analysis (group-ICA) was conducted; a z-score threshold of 2 was applied to these components. The results were acquired in a single group-level ICA using the GIFT toolbox in MATLAB (http://mialab.mrn.org/software/gift/index.html). Voxels in the cerebrospinal fluid, ventricles, eyes, and extracerebral areas (skull) were removed using a brain mask. The ICA was repeated 20 times using ICASSO for stability of the decomposition and determined 13 robust components when combining the 33- and 40 wGA scans, VPT and FT infants combined (Supplementary Figure S6). Among the 13 components, 3 reflected areas related to motion and blood vessels, rather than resting-state activity. These were considered noise and excluded from further analyses. The 10 remaining components were combined into a brain atlas. This group-level ICA-based atlas was manually corrected to ensure anatomical precision. The cerebellum and the brainstem were removed from the final atlas, in order to keep only the cortical and subcortical structures for further analysis, which resulted in 9 RSNs: Limbic (amygdala, insula, anterior cingulate cortex [ACC]), Paralimbic belt (orbitofrontal cortex [OFC], temporal pole), Thalamus, Prefrontal cortex [PFC], Visual, Auditory, Posterior cingulate cortex [PCC], Sensorimotor and Precuneus (Figure 5).

Regional cortical BOLD variability estimation
To extract region-wise time-series from the RS-fMRI data at each time-point, we used our ICA-based brain atlas consisting of 9 RSNs. First, the atlas was registered to each subject’s native time-point specific space (structural T2-weighted) using the diffeomorphic symmetric image normalization algorithm with cross-correlation as similarity metric (SyN-CC), from the Advanced Normalization Tools (ANTs) (Avants et al., 2011) toolbox.
Each subject’s structural T2-weighted image was segmented using the tissue segmentation maps from the UNC neonatal AAL atlas (Shi et al., 2011), which include Gray Matter (GM, cortical and deep gray matter), White Matter (WM) and Cerebral Spinal Fluid (CSF). These maps were first registered to the subject’s native space using the same ANTS algorithm described above, to obtain the three tissue-probability maps (TPMs) in the native space. Next, using these TPMs, the subject’s T2-weighted image was segmented using SPM12, so that subject-specific GM, WM and CSF probability maps were obtained in the subject’s T2 structural space.
To ensure that the BOLD signal was extracted only from GM voxels, the ICA-based brain atlas in the native space was GM-masked using the subject-specific TPMs obtained from the segmentation step. If the probability of one voxel being GM was less than the probability of this voxel being WM or CSF (P(GM) < P(WM) or P(GM) < P(CSF)), then this voxel was not considered as GM and was excluded from the atlas.
The subject’s GM-masked ICA-based atlas in the subject’s structural space was then registered to the subjects’ functional (BOLD) space and the average regional BOLD time-courses were extracted, after performing voxel-wise nuisance regression (CSF, WM, motion signals), detrending and smoothing (Gaussian FWHM 5mm), yielding a matrix with dimensions [#volumes, 9] for each subject and each RS-fMRI time-point (33 and 40 wGA).
The 9 regional time-courses were bandpass filtered ([0.01Hz – 0.1Hz]) to discard noise components and non-resting state fMRI components. BOLD signal variability (BOLD SD) was estimated by computing the sample standard deviation of the time-course of each of the 9 RSNs. This step resulted in 9 regional BOLD variability values per subject and per time-point.
In addition, ALFF (Amplitude of Low-Frequency Fluctuations) was computed by transforming each regional time series into the frequency domain using the Fast Fourier Transform to obtain the power spectrum (Feng et al., 2024; Huang et al., 2020; Yan et al., 2021). The square root of the power spectrum within the 0.01–0.08 Hz (low) frequency band was averaged to obtain regional ALFF values. Regional fALFF (fractional ALFF) was then calculated as the ratio of the low-frequency (0.01–0.08 Hz) ALFF to the total ALFF across all frequencies of the power spectrum. These complementary measures were used to confirm that the results of BOLD SD analyses of spontaneous BOLD fluctuations were consistent across other different metrics estimating BOLD variability.
We selected BOLD variability (BOLD SD) as the primary measure because it captures moment-to-moment fluctuations across the full frequency spectrum, providing a broad and sensitive index of spontaneous neural dynamics. In contrast, ALFF and fALFF are restricted to a predefined low-frequency band and are more dependent on preprocessing choices, whereas BOLD variability provides a more general, assumption-free characterization of spontaneous activity, which is better suited to capturing physiological and microstructural changes during late gestation. ALFF and fALFF results are presented in Supplementary Material.

Multi-Shell Diffusion Imaging pre-processing
MSDI data were preprocessed using MRtrix3 (version 3.0rc3, https://www.mrtrix.org) (Tournier et al., 2019) for denoising, bias field corrections and intensity normalization. The FSL diffusion toolbox (v5.0.11, https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/) (Behrens et al., 2003, Smith et al., 2004) was used for preprocessing of MSDI data. FSL's TOPUP (Andersson et al., 2003, Smith et al., 2004) was used to estimate the off-resonance field, which was then used as input for FSL's EDDY function optimized for neonatal diffusion data, correcting for distortions induced by susceptibility and eddy currents, as well as by motion-induced signal dropout and intra-volume subject movement (Andersson et al., 2016, Andersson et al., 2017, Bastiani et al., 2019). Data were visually inspected to assure quality of motion artifacts correction.

Diffusion cortical microstructural maps
After pre-processing, MSDI data was analysed using three different methodological approaches, each applied to the same set of b-values (b=200, b=1000, b=2000 s/mm2).
To examine cortical microscopic features unconfounded by the directional structure, we use the Spherical Mean Technique (SMT) (Kaden et al., 2016), which focus on the direction-averaged diffusion weighted signal, independent of the fiber orientation distribution. MSDI was analysed using the multi-compartment microscopic diffusion model from the publicly available SMT code (https://github.com/ekaden/smt). SMT metrics corresponding to intra-neurite volume fraction (intra), intrinsic diffusivity (Diff), extra-neurite mean diffusivity (extraMD) and extra-neurite transverse diffusivity (extraTrans) were calculated.
Alongside, we calculated Diffusion Tensor Imaging (DTI) metrics (fractional anisotropy, FA; mean diffusivity, MD) measuring Gaussian diffusion properties of the water in the tissue, as well as Diffusion Kurtosis Imaging (DKI) metrics (mean/radial/axial kurtosis, MK/RK/AK), evaluating Non-Gaussian diffusion and allowing the characterization of tissue microstructural heterogeneity (Jensen & Helpern, 2010; Jensen et al., 2005). For both DTI and DKI metrics, we used the code publicly available from DIPY library (https://docs.dipy.org/stable/interfaces/reconstruction_flow.html), applying the DKI model to the full MSDI data.
The GM-masked ICA-based brain atlas in the subject’s structural space was registered to the subjects’ dMRI space using ANTs (Avants et al., 2011), and averaged diffusion metrics were extracted per cortical RSN and per subject, at each time-point.

Gene expression analysis
Regional patterns of gene expression in the fetal cortex across gestation were obtained from the BrainSpan developmental transcriptomic dataset (http://development.psychencode.org/), comprising mRNA-seq collected from n = 41 specimens aged between 8 pcw and 40 postnatal years. For our study we have included data sampled from 18 post-mortem prenatal brain specimens, aged 8-37 pcw (n = 165 regional tissue samples, mean [SD] age = 17.74 [7.54] pcw, 53% male, mean [SD] postmortem interval = 6.52 [11.03] hours, mean [SD] RNA integrity number [RIN] = 9.34 [0.66]).
Tissue had been collected after obtaining parental or next of kin consent and with approval by the institutional review boards at the Yale University School of Medicine, the National Institutes of Health, and at each institution from which tissue specimens were obtained.
For each brain, regional dissection including 11 neocortical regions was performed (dorsolateral PFC [DLPFC], ventrolateral PFC [VLPFC], OFC, medial frontal cortex [MFC], primary motor/sensory/auditory/visual cortex [M1, S1, A1C, V1], inferior parietal cortex [IPC], superior temporal cortex [STC] and inferior temporal cortex [ITC]). Tissue processing and detailed anatomical boundaries for each cortical region at each stage of development are provided elsewhere (Kang et al., 2011; Li et al., 2018). Regional tissue samples were subject to mRNA-seq using an Illumina Genome Analyzer IIx (Illumina, San Diego, California, United States of America) and mRNA-seq data processed using RSEQtools (version 0.5) (Habegger et al., 2011). Gene expression was measured as RPKM. Conditional quantile normalization was performed to remove guanine-cytosine (GC) content bias, and ComBat was used to remove technical variance due to processing site (Yale or University of Southern California) (Hansen et al., 2012; Johnson et al., 2007; Li et al., 2018).
The prenatal gene expression data were initially filtered to only include protein-coding genes (NCBI GRCh38.p12, n = 18,524 out of a possible 20,720). In order to restrict our analysis to focus on genes expressed in the developing cortex, this list was further filtered to only contain genes expressed by cells in the fetal cortex based on the composite list of prenatal cell markers from 5 independent single-cell RNA studies of the developing fetal cortex (Ball et al., 2020). This resulted in expression data from a final set of 5,287 genes.
For our study, gestational periods were grouped into early-fetal (T0, 8-13 pcw), mid-fetal (T1, 16-22 pcw) and late-fetal (T2, 35-37 pcw). The neocortical regions with available genetic data were matched to our RSNs as follows: MFC = Limbic (1), OFC = Paralimbic belt (2), DLPFC+VLPFC = PFC (4), V1 = Visual (5), A1C = Auditory (6), M1+S1 = Sensorimotor (8), IPC = Precuneus (9). Next, these regions were grouped according to the presence (group 1) or absence (group 2) of significant BOLD variability changes from 33- to 40-wGA (Figure 3a-b).

Neonatal and demographic data statistical analysis
To test for differences in demographic and perinatal data between VPT and FT groups, categorical variables (sex, intrauterine growth restriction, sepsis, intraventricular hemorrhage grade I) were analysed using chi-squared test, whereas continuous variables were compared using independent samples t-test, with group as an independent variable and the following dependent variables: GA at birth, GA at MRI at TEA, birth weight, birth height, birth head circumference, APGAR score at 1 and 5 min after birth and parental socio-economic status score (SES) (Largo et al., 1989).

MRI data statistical analysis
Longitudinal (from 33- to 40-wGA) and cross-sectional differences (at TEA, between VPT at 40-wGA and FT newborns) in both cortical BOLD variability and microstructural diffusivities were assessed using two-sided paired-samples t-tests and independent samples t-tests, respectively, with FDR correction for multiple comparisons. Prior to statistical testing, outliers were identified for each metric, by region and timepoint, based on 1.5×interquartile range (IQR). Values falling below Q1−1.5×IQR or above Q3+1.5×IQR were excluded from subsequent analyses. All analyses and visualizations were performed using RStudio (R version 4.4.1).
For visualization of inter-individual variability using boxplots, all subject-level data were normalized by the maximum absolute change across all participants and regions for each metric, thereby constraining values to the range [−1,1] while preserving the directionality of effects. These globally normalized values were used to generate boxplots, enabling comparison of relative regional patterns across participants and across metrics, independent of absolute effect magnitude. Additional boxplots displaying the raw, non-normalized data, are provided in Supplementary Material (Figure S2).
For bar-plot visualizations and the clustering analysis, for each metric, mean regional changes were first computed for each cortical resting-state network region and then normalized per metric by dividing by the largest absolute mean regional change, preserving the directionality of effects. This scaling yielded a group-level region × metric matrix representing the relative spatial distribution of longitudinal changes across cortical resting-state network regions.
For Lasso-based regression modelling, the same normalized group-level region × metric matrix was used (mean regional changes were first computed for each cortical resting-state network region and then normalized per metric by dividing by the largest absolute mean regional change), but expressed as absolute values, in order to focus on the magnitude of change, irrespective of the direction. This approach prevents potential cancellation of positive and negative effects across regions and ensures that predictors with strong effects in opposite directions are appropriately retained.
Microstructural diffusivities were clustered using the R-package and function ConsensusClusterPlus (k-means base algorithm and Euclidean distance) (Wilkerson & Hayes, 2010), to identify stable clusters of regional diffusivities. The optimal number of clusters was determined using the elbow method, based on the within-cluster sum of squares (WCSS). This approach enabled grouping cortical regions based on the similarity of their longitudinal microstructural maturation, and to compare these microstructural clusters with those derived from BOLD variability longitudinal changes.
Lasso-based sparse multivariate regression models were used to identify the most relevant microstructural predictors of BOLD variability longitudinal changes, between 33- to 40-wGA. Analyses were performed using the glmnet package in R with leave-one-region-out cross-validation. The analysis was conducted on the group average 9 region x 10 metric matrix (n = 9 regions), where BOLD variability change served as the dependent variable and nine microstructural metrics (extraMD, extrTrans, diff, intra, MD, FSA, MK, AK, RK) served as predictors. All metrics were standardized prior to modelling, as previously described. Model selection was performed using leave-one-region-out cross-validated mean squared error, and for each penalty value, RMSE and R²ₗₒₒcᵥ were calculated. Two penalty choices were evaluated: the more conservative penalty corresponding to one standard error above the minimum (λ₁se) and the less penalized minimum-error solution (λ_min). In addition to a pure LASSO model (α = 1), a LASSO Elastic Net model (α = 0.5) was examined to allow for partial retention of correlated predictors. Predictors with non-zero coefficients at the selected penalty values were considered relevant.

Gene expression data analysis
A linear mixed effects regression model (LMM) was used to identify significant changes in gene expression as a function of regional group (1 and 2) and time, with fixed effects of sex and RIN (RNA integrity number), including specimen ID as a random effect. We identified genes with a significant group *time interaction and calculated estimated marginal means, followed by pairwise comparisons to assess differences in gene expression over time and between groups. Multiple comparisons were adjusted using false discovery rate (FDR) correction. Modelling was performed in RStudio (R version 4.4.1) using “lmer()” function from lme4 package. Gene function was classified using gene set enrichment analysis (GSEA) using WebGestalt (Wang et al., 2017). Enrichment ratios were calculated as the proportion of cell class-specific genes in our gene list of interest compared to the proportion in the full background set. The background gene set was defined as the full list of protein-coding genes included in the analysis (n = 5,287 expressed by cells in the fetal cortex). We corrected for multiple comparisons across cell classes using FDR. The genes with a significant group*time interaction, found to have a significant enrichment ratio linked to neocortical organization, were combined for expression analysis across time points (T0, T1 and T2) and between groups (1 and 2) using two-way ANOVA followed by post-hoc comparisons using Tukey's test (HSD).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The non-invasive imaging readout - BOLD signal variability - indirectly reflects spontaneous neuronal activity and varies across cortical regions, cognitive tasks, and age. This manuscript aims to investigate how BOLD variability evolves in very preterm (VPT) infants, how it correlates to cortical microstructure and gene expression, and how it differs from full-term (FT) newborns at term-equivalent age (TEA).\nTo achieve this, the authors collected resting-state fMRI and multi-shell diffusion imaging for VPT infants (at 33 and 40 weeks GA) and FT infants (40 weeks GA). They first compared BOLD variability longitudinal changes in VPT infants between 33 and 40 weeks GA. Primary sensory, sensorimotor and pDMN regions (cluster 1) had increased BOLD variability across 33 to 40 weeks. The authors then correlated longitudinal changes in BOLD variability with cortical microstructural changes as measured by diffusion imaging metrics. In regions where BOLD variability increased with GA, diffusivity decreased – revealing a relationship between cortical microstructure and BOLD variability. Next, the authors associated longitudinal microstructural and BOLD variability changes with gene expression changes using the brainspan developmental transcriptome. This analysis revealed that cluster 1 regions had increased gene expression profiles associated with neuronal ensheathment and gliogenesis. Lastly, they compared BOLD variability and cortical microstructure at term equivalent age between VPT and FT infants revealing regions-specific decreases in BOLD variability and global increases in diffusivities in VPT infants at TEA.\nOverall this study uses state-of-the-art methodologies to uncover the biological significance of BOLD variability changes in the early development of VPT infants. This is an important topic as the field would benefit from a better understanding of underlying mechanisms of preterm birth-induced brain injuries. However, the study has several issues that should be investigated and resolved.\n1)The applicability of fetal gene expression patterns associated with cluster 1 regions to the longitudinal changes in BOLD variability is difficult to interpret. Surely, there is an increase in gene expression associated with gliogenesis and myelination from T1 to T2, however, this is a much broader developmental timeframe (16-22GA to 35-37GA) than is assessed with imaging (33GA to 40GA).\n2)Gene expression was analyzed using a filtered list of protein-coding genes that were identified in single-cell datasets of the fetal cortex. Single-cell transcriptomic techniques often have low sensitivity and do not fully capture RNAs expressed at low levels. The filtering performed here greatly limits the potential to identify novel biological processes that may be underlying physiological changes identified here.\n3)The novelty and overall significance of these findings are difficult to assess without contextual discussion of two recent studies examining resting-state fMRI signals in preterm infants (Drayne et al., 2024 & Mella et al., 2024).\n4)TEA and FT imaging data is taken from infants that have been exposed to postnatal environments for ~7 weeks vs. a few days, respectively. The analysis of BOLD variability and microstructural differences between these groups do not consider the role of exogenous factors that may influence their interpretation of the data.", "reviewer_2": "This manuscript aims to characterize changes in BOLD variability and their structural and biological underpinnings in very preterm (VPT) infants, with comparisons to full-term (FT) infants. The authors used fMRI to assess regional BOLD variability, diffusion MRI to evaluate cortical microstructure across longitudinal preterm data (33w GA to TEA), and independent gene expression analyses to explore biological correlates of BOLD variability. A LASSO model was applied to investigate the relationship between microstructure and functional variability. It is very interesting to study the difference between VPT and FT from the perspective of BOLD signal variability and its relationship with cortical microstructure and gene expression. Addressing the interplay between BOLD variability from resting-state fMRI, microstructure measures from diffusion MRI, and gene expression is a major strength of this work. However, there are notable concerns about the fundamental mechanisms of BOLD variability and rigorousness of methodology, which limit the impact and robustness of the current work.\nMajor Comments\n1. The fundamental mechanisms underlying BOLD variability are not clear. Thus, it does not appear robust to establish a study centered on a measure with its mechanism not well understood. The limited citations provided to explain BOLD variability are constrained to few groups, raising the question if BOLD variability is broadly reproducible or widely applicable.\n2. Methodology needs to be further improved to be more rigorous. For example, the description of the LASSO model is relatively vague. It is not clear how the regression was performed, how significance was determined, or what the model performance was. In addition to listing the selected three measures, the manuscript should report the actual sample size, model fit statistics, and p-values for all predictors. As another example, the analysis of cortical gray matter microstructure using diffusion MRI should consider mitigating partial volume effects.\n3. The relatively small sample size, combined with the lack of validation on publicly available datasets, raises concerns about the generalizability of the findings. The results would be more convincing with rigorous validation. Although the authors note the difficulty of acquiring longitudinal preterm participants, the Developing Human Connectome Project (dHCP) contains a large cohort of both preterm and term infants, including longitudinal scans. Leveraging this public dataset would greatly strengthen the conclusions.\n4. While the authors chose 9 microstructure measures from diffusion MRI, BOLD variability was used as the only measure derived from resting state fMRI. Other relatively more widely adopted measures from resting-state fMRI such as amplitude of low frequency fluctuations (ALFF) or regional homogeneity (ReHo) should be included to test the findings or the authors should justify exclusion of these other more widely used measures from resting state fMRI. Lacking these measures from resting state fMRI weakens the functional analysis and limits the ability to fully address how microstructural and biological development support functional maturation.\nMinor Comments\n1. While the Introduction provides useful background on diffusion models, it does not sufficiently justify why DTI measures (e.g., FA, MD) were applied to study cortical microstructure, given their traditional use in white matter.\n2. Since the authors used three different diffusion models, it would be beneficial to specify the b-values used in each diffusion model in the Methods.\n3. The definition of “clusters” is not clear. In BOLD variability, “clusters” appear to denote groups of regions showing significant increases between two timepoints. But in diffusion MRI, k-means clustering (k=3) was used without justification of why three clusters were chosen. To avoid confusion, the authors could consider using “groups” for the BOLD variability groups of regions and provide validation (e.g., elbow method, silhouette score) for the k-means choice of clusters.\n4. The authors might want to be cautious about using “trajectory”, as only two timepoints were included in the manuscript.\n5. There is a typo in the last sentence of the Introduction: “weather” should be corrected to “whether.”", "reviewer_3": "Thank you for asking me to review this interesting paper examining neonatal brain development.\nThe work combines MRI data with open access transcriptomic data to produce interesting and plausible biologic insights. The combination of BOLD as a measure of cortical function with transcription data is as far as I know completely novel and of great interest. The combination of MRI with transcription data is an emerging and promising area of science and Ball is the world leader in developing this approach in the developing brain. The paper leverages longitudinal imaging data encoding the 'experiment of nature' that is preterm delivery. Premature birth is a wide-ranging and robust disrupter of brain development, which has been shown (by these authors as well as many others) to affect almost every aspect of brain development and so provides a useful 'positive control' reality check on these results.\nThe imaging data appropriately collected with reasonably up to date sequences. The cohort is small, and indeed quite a lot smaller than most substantial studies currently being published. The data were collected by the researchers, and they probably trust their own data quality, relying on the longitudinal data to add value to the small cohort. Nevertheless, in the same way that they used open access transcriptomic data, they could have replicated their results using the open access imaging data now freely available, which are an order of magnitude larger than the current dataset, and which some of this group have used in the past.\nThe image analysis is reasonable. It appears that there was no motion correction applied to the scan data, but a reasonable cut-off for exclusion of scans was used to remove corrupted datasets. It is a bit controversial whether this scrubbing of individual data (which disrupts the time course) affects fMRI interpretation, but my view is that it is unlikely to invalidate measures of BOLD variability. Variability is a nice and under-explored measure of cortical function, the results are interesting and as far as I know novel, and the functional locations revealed are interesting.v Spherical means is a good method for analysing diffusion data and the results are consistent with previous studies using different approaches.\nThe imaging results are plausible and interesting, with the increases in visual and sensory-motor cortex in late gestation as expected given known developmental pathways, and to is also interesting that they did not find BOLD variability increases in pre-frontal cortext or limbic- it is plausible to interpret this as the frontal cortex not having started its rapid development while the limbic system may have completed its comparable phase. However it would have been interesting to see more data on temporal cortex which is rapidly developing at this time.\nThe transcriptomic data and analysis are high quality and appropriate. it is helpful to have the three time points in the data even if there are no imaging data to compare to the first period. The gene set analysis is uncontroversial and Lasso is a reasonable approach to this underdetermined dataset. As with the imaging data the results are plausible and it is re-assuring that neuronal ensheathment is most significant relation, which gives confidence in other pathways detected. It is a shame that the analysis does not make use of the remarkable microbrain histologic parcellation that one of the authors has recently published, but relies on small patch location for the transcriptomic analysis, but that does not invalidate the results. The gene sets exposed in the analysis are interesting and add to existing knowledge."}, "author_response": "We thank sincerely the reviewers and the editor for their constructive comments and suggestions that have greatly improved our manuscript. Detailed responses are provided below in blue, and additions to the manuscript are highlighted in yellow.\nReviewer #1:\nOverall, this study uses state-of-the-art methodologies to uncover the biological significance of BOLD variability changes in the early development of VPT infants. This is an important topic as the field would benefit from a better understanding of underlying mechanisms of preterm birth-induced brain injuries. However, the study has several issues that should be investigated and resolved.\n1. The applicability of fetal gene expression patterns associated with cluster 1 regions to the longitudinal changes in BOLD variability is difficult to interpret. Surely, there is an increase in gene expression associated with gliogenesis and myelination from T1 to T2, however, this is a much broader developmental timeframe (16-22GA to 35-37GA) than is assessed with imaging (33GA to 40GA).\nWe agree with the reviewer that the developmental time frame is not exactly the same between the fetal gene expression and our neuroimaging data. The available developmental transcriptomic dataset (http://development.psychencode.org/) is limited by the inherent difficulty in obtaining well-characterized human fetal brain samples across the full gestational period. It provides data sampled from 18 post-mortem prenatal brain specimens aged 8-37 pcw, which we have clustered as proposed in (Kang et al., 2011), according to developmental windows. Since only two specimens fall within the specific developmental window in which our cohort was scanned (Figure R1A), excluding samples from the broader prenatal period would prevent us from modelling the temporal gene-expression trends occurring during our observational window. We would also like to clarify that, as we expect gene expression changes to precede structural and functional changes, the available timeframes align well with the aims of our study, allowing us to model change in gene expression that occur both before and during the observed regional differences in functional and structural maturation in the cortex.\nOur results are consistent with the findings from (Kang et al., 2011), illustrating the well-described temporal sequence in cortical gene expression, with cell proliferation genes peaking first, followed by expression of immature neuron markers (e.g., DCX, marker of neurogenesis and neuronal migration), and subsequently gliogenesis and myelination gene programs (Figure R1C, (Kang et al., 2011)). The expression of genes involved in cortical myelination increases well before myelin becomes detectable on histology, which is typically after birth (Figure R1B, (Miller et al., 2012)). In line with Kang et al., we find a marked increase in myelination-related gene expression beginning from the mid-fetal period (100 days post-conception ≈ 16wGA) and continuing until term age. We show that from 33 to 40 wGA, while myelination-related gene expression is rising, regions exhibiting the largest increases in BOLD variability in our dataset correspond to those with a significantly higher expression in those genes. Plus, the expression of these genes remains significantly higher in these same regions during the late-fetal period, compared to regions without significant BOLD variability changes. These findings support the idea that early transcriptional programs lay the groundwork for later functional changes.\nWe have added this rational into the Discussion, section “Genes mediating gliogenesis and neuronal ensheathment increase more in RSNs where BOLD variability increased”, as follows:\nBetween 16-22 pcw and 35-37 pcw, corresponding to the mid- to late-fetal period, analysis of ex-vivo fetal cortical gene expression patterns revealed a significant upregulation of genes associated with neuronal ensheathment and gliogenesis. These findings are consistent with previous literature showing a marked increase in myelination-related gene expression starting from the mid-fetal period, continuing until birth, and into post-natal period (Kang et al., 2011).\nThis increase in genes associated with neuronal ensheathment and gliogenesis, from mid- to late-fetal period, is significantly higher in the regions where BOLD variability increases from 33- to 40-wGA, namely the primary sensory (visual and auditory), sensorimotor and proto-DMN (precuneus and PCC). The expression of these genes remains significantly higher in these same regions during the late-fetal period, compared to regions without significant BOLD variability changes. Our findings support the idea that early transcriptional programs precede and lay the groundwork for later functional changes in the developing brain.\nFigure R1. A. Bar plot showing the number of observations (specimen ID) per gestational age (in weeks), obtained from the BrainSpan developmental transcriptomic dataset (http://development.psychencode.org/). B. Adapted from (Miller et al., 2012), bar graph depicts mean percent of maximum mature adult myelinated fiber length density (MFLD) across development in humans, across cortical regions. C(d). Adapted from (Kang et al., 2011), sets of genes (expressed as percentage of maximum) plotted against age to represent general trends and regional differences in several neurodevelopmental processes in NCX (neocortex), HIP (hippocampus) and CBC (cerebellum).\n2. Gene expression was analyzed using a filtered list of protein-coding genes that were identified in single-cell datasets of the fetal cortex. Single-cell transcriptomic techniques often have low sensitivity and do not fully capture RNAs expressed at low levels. The filtering performed here greatly limits the potential to identify novel biological processes that may be underlying physiological changes identified here.\nWe thank the reviewer for their comment. We agree that single-cell transcriptomic datasets may have reduced sensitivity for low-abundance transcripts. We emphasize that BrainSpan (http://development.psychencode.org/) which provides bulk mRNA-seq data, was used as the primary dataset, and that the single-cell datasets were used solely to define a high-confidence list of genes expressed in the developing cortex, not as the source of the expression profiles analyzed in our study.\nThe BrainSpan bulk mRNA-seq dataset was generated using standardized protocols with stringent quality-control criteria, aiming to minimize false positives and capture robust region-level transcriptional patterns. Our filtering strategy follows standard practice for restricting bulk transcriptomic analyses to genes with validated prenatal cortical expression, using high-confidence markers derived from multiple independent single-cell datasets, to ensure that our analysis targeted biologically meaningful cortical signals. This consensus reduces noise from transcripts not expressed in the fetal cortex, increases specificity for robust developmental processes, and minimizes false positives introduced by bulk-level technical variation. Our aim was to characterize developmental patterns in genes with validated relevance to fetal cortical biology, rather than to conduct an exhaustive exploratory screen for novel low-abundance transcripts.\nThe developmental trajectories identified with this biologically grounded gene set are consistent with previous literature (Miller et al., 2012) and, importantly, align with the functional and structural neuroimaging changes observed during this period, indicating that the major biological processes are robustly captured despite this conservative filtering.\n3. The novelty and overall significance of these findings are difficult to assess without contextual discussion of two recent studies examining resting-state fMRI signals in preterm infants (Drayne et al., 2024 & Mella et al., 2024).\nWe thank the reviewer for highlighting these interesting recent studies. Indeed, Drayne et al. (2024) and Mella et al. (2024) examine resting-state fMRI signals in preterm infants using the Hurst Exponent (HE) to characterize the temporal complexity of BOLD signals. While both HE and BOLD variability (BOLD SD) capture aspects of neural signal dynamics, the two metrics reflect distinct properties. BOLD variability (BOLD SD) specifically quantifies the moment-to-moment fluctuations in signal amplitude, providing a measure of signal variability sensitive to region-specific functional activity. In contrast, HE captures temporal correlations over longer time scales, aiming to characterize the temporal complexity of BOLD signals, indexing the persistence or “memory” (how past fluctuations of the signal predict future fluctuations) of the signal, rather than its instantaneous amplitude variability.\nAs suggested, we have added to our discussion the results from Drayne et al. (2024) and Mella et al. (2024) studies, adding these independent but related perspectives on functional maturation in preterm infants.\nDiscussion, Section “Cortical microstructural maturation underlies RSNs BOLD variability increases during preterm infants’ brain development”:\nOur findings, revealing significant longitudinal increases in BOLD variability specifically in primary motor, sensory, and posterior brain regions from 33wGA to TEA, with no changes in prefrontal regions, agree with the spatiotemporal pattern of synaptogenesis and cortical activity. These results, quantifying fluctuations in BOLD signal amplitude, complement previous literature findings showing that the temporal complexity of BOLD signals, evaluated by means of the Hurst Exponent (HE), also increases overall in very preterm infants’ gray-matter from birth until TEA (Drayne et al., 2024; Mella et al., 2024). Similar to our findings, the greatest increases in HE were also found in motor and visual networks, compared to frontal networks, suggesting that primary sensorimotor areas mature from “more random” to “more ordered” states faster than the frontal ones (Drayne et al., 2024). BOLD variability provides thus a complementary perspective, suggesting that primary sensorimotor regions also present increasingly variable functional activity, highlighting their accelerated functional maturation.\nDiscussion, Section “Preterm birth impacts BOLD variability and cortical microstructure at TEA”:\nCompared to FT newborns, very premature infants at TEA showed reduced BOLD-variability specifically in sensory RSNs (auditory and visual), PCC, thalamus and limbic, denoting the impact of preterm birth on cortical functional maturation. Our findings further complement recent evidence that the Hurst Exponent, which measures the temporal complexity of BOLD signals, is similarly reduced in somatosensory-motor, visual and auditory networks in preterm infants at TEA, compared to FT (Mella et al., 2024).\n4. TEA and FT imaging data is taken from infants that have been exposed to postnatal environments for ~7 weeks vs. a few days, respectively. The analysis of BOLD variability and microstructural differences between these groups do not consider the role of exogenous factors that may influence their interpretation of the data.\nWe thank the reviewer for raising this point. The primary aim of our analysis was to characterize whether VPT infants at TEA (and thus exposed to a postnatal environment early) differ from FT infants in cortical BOLD variability and microstructure. Because our goal was to characterize group-level differences rather than infer causal relationships, we applied inferential statistics without covariates. Importantly, early postnatal environmental exposure is an intrinsic feature of preterm birth. Adjusting this analysis for additional specific clinical exposures of the VPT group might inadvertently remove biologically meaningful variance that is not necessarily confounding to the contrast of interest.\nWe would also like to highlight that our VPT cohort was relatively homogeneous: included only infants born before 32 weeks of gestation who were clinically stable and free of major complications. Notably, there were no statistically significant differences between the VPT and FT groups in demographic and clinical variables such as sex, socioeconomic status (SES), intrauterine growth restriction (IUGR), bronchopulmonary dysplasia (BPD), or intraventricular hemorrhage (IVH). The GA at the MRI scan was also similar (Table 1). This strengthens the interpretability of the observed group differences.\nWhile the range of potential exogenous factors possibly playing a role during this developmental period will likely be of interest to researchers in the field, we feel it is beyond the scope of this study to explore these factors in detail, deferring to future research using an experimental protocol designed specifically to do so.\nFinally, our findings show that although VPT infants exhibit increasing BOLD variability and progressive cortical microstructural maturation from preterm birth to TEA, by TEA they nonetheless demonstrate reduced functional and microstructural cortical maturation compared with FT newborns, reinforcing VPT infants’ distinct developmental trajectory.\nReviewer #2:\nIt is very interesting to study the difference between VPT and FT from the perspective of BOLD signal variability and its relationship with cortical microstructure and gene expression. Addressing the interplay between BOLD variability from resting-state fMRI, microstructure measures from diffusion MRI, and gene expression is a major strength of this work. However, there are notable concerns about the fundamental mechanisms of BOLD variability and rigorousness of methodology, which limit the impact and robustness of the current work.\nMajor Comments\n1. The fundamental mechanisms underlying BOLD variability are not clear. Thus, it does not appear robust to establish a study centered on a measure with its mechanism not well understood. The limited citations provided to explain BOLD variability are constrained to few groups, raising the question if BOLD variability is broadly reproducible or widely applicable.\nWe thank the reviewer for their comment, and we agree that the fundamental mechanisms underlying BOLD variability are not yet completely understood. However, converging evidence from multiple independent research groups supports its interpretive value as a proxy for spontaneous neural dynamics and network-level organization. We have revised the introduction to provide additional background and include further relevant references (changes, in yellow, can be seen below and in the manuscript).\nSeveral independent groups have investigated BOLD variability across different populations and age ranges using complementary methodological approaches. In fact, BOLD variability has been operationalized using slightly different metrics, including standard deviation (SD), mean squared successive difference (MSSD), as well as the closely related amplitude of low-frequency fluctuations (ALFF/fALFF). Although ALFF focuses on a specific low-frequency band (0.01–0.08 Hz), whereas SD and MSSD summarize overall temporal variance, the dominant power of the BOLD signal lies within the low-frequency range (0.01-0.1Hz). All these metrics are mathematically related and empirically correlated. Zhang et al. (2022) and Nomi et al. (2017) demonstrated that SD and MSSD are highly correlated in resting-state fMRI data, and ALFF and SD approaches have been recognized as capturing similar aspects of BOLD signal fluctuations over time (Lim et al., 2021; Millar et al., 2020). Therefore, these metrics capture very similar spontaneous BOLD fluctuations, which reflect underlying neural dynamics.\nFindings from multiple research groups demonstrate that BOLD variability is a widely applicable and reproducible measure.\nWe present a table below summarizing the metrics used by different groups and the populations studied. The limited number of studies particularly in developmental populations, highlights the need for more research during this period. Our study further extends its application to very preterm infants before term-equivalent age (prior to the typical time of birth), a population that has not been studied using these metrics. We have included additional references from various groups in the introduction and discussion, when relevant, as shown below and in the manuscript.\nIntroduction\nHowever, the “BOLD signal variability,” characterized by its variance, captures the magnitude of moment-to-moment regional variations in BOLD signal (Garrett et al., 2010, 2011, 2013). It is thought to indirectly reflect spontaneous brain activity arising from moment-to-moment fluctuations in neuronal activity. In fact, BOLD variability has been shown to correlate with electrophysiological measures of neural dynamics, including EEG power and temporal signal variability (Baracchini et al., 2021; McIntosh et al., 2010; Misic et al., 2011; Uddin, 2020; Vakorin et al., 2011). These temporal fluctuations in neuronal activity, as measured by electrophysiological and EEG studies, are thought to contribute to synaptic connectivity and associate with better cognitive performance (Faisal et al., 2008; McIntosh et al., 2008; Misic et al., 2011; Uddin, 2020; Vakorin et al., 2011), thereby underlying brain function.\nAlthough BOLD variability is affected by vascular factors such as baseline cerebral blood flow and neurovascular coupling (Whittaker et al., 2016), studies combining BOLD and arterial spin labeling (ASL) imaging show that group differences remain even after accounting for these vascular influences (Garrett et al., 2017; Roberts R.P., 2023). This evidence indicates that BOLD variability reflects not only vascular effects but also neural dynamics, supporting its use as a meaningful index of underlying functional brain activity.\nIn adults, BOLD variability varies across cortical regions, with greater variability in association and transmodal regions (Baracchini et al., 2021). Greater BOLD variability has been linked to increased network organization, greater cognitive performance and increased task difficulty across different studies (Baracchini et al., 2021; Baracchini et al., 2023; Garrett et al., 2011, 2013; Goodman et al., 2024; Nomi et al., 2017; Protzner et al., 2013; Roberts R.P., 2023; Zhang et al., 2022), highlighting its functional relevance.\nBOLD variability typically increases during childhood (Wang et al., 2021) and serves as a strong predictor of age across the lifespan (Garrett et al., 2010), then declining with aging (Baracchini et al., 2023; Garrett et al., 2013; Grady & Garrett, 2014; Rieck et al., 2022). To date, no studies have examined changes in BOLD variability in infants prior to term-equivalent age (TEA).\nDiscussion, Section “Preterm birth impacts BOLD variability and cortical microstructure at TEA”:\nCompared to FT newborns, very premature infants at TEA showed reduced BOLD-variability specifically in sensory RSNs (auditory and visual), PCC, thalamus and limbic, denoting the impact of preterm birth on cortical functional maturation. These results align with previous work showing that BOLD variability, measured as ALFF (amplitude of low frequency fluctuations), is reduced in moderate-to-late preterm infants at TEA, compared to FT newborns, predominantly in the primary sensory and motor cortices, and in the posterior cingulate cortex and precuneus (Wu et al., 2016).\nTable R1. Non-exhaustive summary table of metrics used across groups (with full citations at the end):\nMetric Study Group Population\nBOLD SD (Garrett et al., 2010, 2011, Adults & Older adults\n2013; Garrett et al., 2017)\n(Roberts R.P., 2023) Adults\n(Goodman et al., 2024) Lifespan\n(Protzner et al., 2013) Adults\n(Millar et al., 2020) Adults and older adults\n(Lim et al., 2021) Adults\n(Scarapicchia et al., 2018) Older adults\nMSSD (Zhang et al., 2022) Young adults\n(Mean Squared MSSD & BOLD SD\nSuccessive Difference)\n(Nomi et al., 2017) Lifespan\nMSSD & BOLD SD\n(Baracchini et al., 2021; Adults\nBaracchini et al., 2023)\n(Easson & McIntosh, 2019) Children and adolescents\n(Boylan et al., 2021) Adults and older adults\nALFF / fALFF (Wu et al., 2016) Moderate-to-late preterm\ninfants\n(Huang et al., 2020) Neonates vs adults\n7\n(Zang et al., 2007) Children\n(Zou et al., 2008) Children and adults\n(Zhao et al., 2018) Young adults\n(Long et al., 2017) Children\n(Montala-Flaquer et al., 2022) Adults and older adults\n(Lee & Hsieh, 2017) Adults\n(Zheng et al., 2019) Older adults\n2. Methodology needs to be further improved to be more rigorous. A) For example, the description of the LASSO model is relatively vague. It is not clear how the regression was performed, how significance was determined, or what the model performance was. In addition to listing the selected three measures, the manuscript should report the actual sample size, model fit statistics, and p-values for all predictors. B) As another example, the analysis of cortical gray matter microstructure using diffusion MRI should consider mitigating partial volume effects.\nWe thank the reviewer for their comment. We have revised overall the Results section, adding additional details where necessary, as highlighted in the manuscript.\nA)\nIn particular, we have now clarified the LASSO analysis in the Methods section, as described below in yellow and in the manuscript.\nIn resume, we performed the LASSO regression on the group-average 9 region × 10 metric matrix (n = 9 regions) to identify the microstructural measures most strongly associated with regional changes in BOLD variability. This region-level approach was chosen to emphasize biologically interpretable spatial patterns, smoothing out inter-subject variability and noise. Subject-level longitudinal analyses would require accounting for the non-independence of multiple regions within the same subject. Since regions from the same subject are correlated (not independent), performing LASSO at the subject level would violate regression assumptions, which is why we opted for the group-average approach.\nRegarding significance testing, LASSO is a penalized regression and feature-selection method. Therefore, classical p-values for individual coefficients are not defined or appropriate, as penalization biases coefficient estimates and violates standard inference assumptions. Instead, variable importance is determined by whether the predictors retain non-zero coefficients at the cross-validated λ (the LASSO penalty parameter), rather than by hypothesis testing. In this study, λ was selected using cross-validation, and model generalizability was explicitly evaluated using leave-one-region-out cross-validation. Reported regression coefficients therefore reflect relative feature importance, rather than statistically independent effects. This is consistent with the statistical literature and with standard practice.\nWe acknowledge that the Lasso-based results have changed slightly following revision, given the implementation of leave-one-region-out cross-validation, as well as an identified error in outlier selection for some metrics prior to computing regional averages. Importantly, these updates do not alter the main findings or overall interpretation. All analyses are fully reproducible using the “Source Data” provided with the manuscript, along with the R code used for all analyses, which will be made publicly available on Zenodo upon publication.\nMethods, Section “MRI data analysis”:\nLongitudinal (from 33- to 40-wGA) and cross-sectional differences (at TEA, between VPT at 40-wGA and FT newborns) in both cortical BOLD variability and microstructural diffusivities were assessed using paired-samples t-tests and independent samples t-tests, respectively, with FDR correction for multiple comparisons. Prior to statistical testing, outliers were identified for each metric, by region and timepoint, based on 1.5×interquartile range (IQR). Values falling below Q1−1.5×IQR or above Q3+1.5×IQR were excluded from subsequent analyses. All analyses and visualizations were performed using RStudio (R version 4.4.1).\nFor visualization of inter-individual variability using boxplots, all subject-level data were normalized by the maximum absolute change across all participants and regions for each metric, thereby constraining values to the range [−1,1] while preserving the directionality of effects. These globally normalized values were used to generate boxplots, enabling comparison of relative regional patterns across participants and across metrics, independent of absolute effect magnitude. Additional boxplots displaying the raw, non-normalized data, are provided in Supplementary Material (Figure S2).\nFor bar-plot visualizations and the clustering analysis, for each metric, mean regional changes were first computed for each cortical resting-state network region and then normalized per metric by dividing by the largest absolute mean regional change, preserving the directionality of effects. This scaling yielded a group-level region × metric matrix representing the relative spatial distribution of longitudinal changes across cortical resting-state network regions.\nFor Lasso-based regression modelling, the same group-level region × metric matrix was used, but expressed as absolute values, in order to focus on the magnitude of change, irrespective of the direction. This approach prevents potential cancellation of positive and negative effects across regions and ensures that predictors with strong effects in opposite directions are appropriately retained.\nMicrostructural diffusivities were clustered using the R-package and function ConsensusClusterPlus (k-means base algorithm and Euclidean distance) (Wilkerson & Hayes, 2010), to identify stable clusters of regional diffusivities. The optimal number of clusters was determined using the elbow method, based on the within-cluster sum of squares (WCSS). This approach enabled grouping cortical regions based on the similarity of their longitudinal microstructural maturation, and to compare these microstructural clusters with those derived from BOLD variability longitudinal changes.\nLasso-based sparse multivariate regression models were used to identify the most relevant microstructural predictors of BOLD variability longitudinal changes, between 33- to 40-wGA. Analyses were performed using the glmnet package in R with leave-one-region-out cross-validation. The analysis was conducted on the group average 9 region x 10 metric matrix (n = 9 regions), where BOLD variability change served as the dependent variable and nine microstructural metrics (extraMD, extrTrans, diff, intra, MD, FA, MK, AK, RK) served as predictors. All metrics were standardized prior to modelling, as previously described. Model selection was performed using leave-one-region-out cross-validated mean squared error, and for each penalty value, RMSE and R²ₗₒₒcᵥ were calculated. Two penalty choices were evaluated: the more conservative penalty corresponding to one standard error above the minimum (λ₁ se) and the less penalized minimum-error solution (λ_min). In addition to a pure LASSO model (α = 1), a LASSO Elastic Net model (α = 0.5) was examined to allow for partial retention of correlated predictors. Predictors with non-zero coefficients at the selected penalty values were considered relevant.\nResults, Section “Relationship between longitudinal cortical microstructural changes and BOLD variability delta changes”:\nTo identify the microstructural metrics that best predicted the observed changes in BOLD variability, we employed Lasso-based sparse multivariate regression models on the regional-average data. Out of the 9 microstructural metrics, 4 features (MD, RK, AK and FA) were retained, at the optimal pure LASSO penalty (λ₁ se = 0.015, α = 1), as predictors of regional BOLD variability changes from 33- to 40-wGA. MD showed the strongest association (β = 0.419), followed by RK (β = 0.402), AK (β = 0.353), and FA (β = 0.332). The model achieved a high training fit (R² = 0.965); however, leave-one-region-out cross-validation indicated that these longitudinal changes explained 28% of the observed variance in BOLD variability (R²ₗₒₒcᵥ = 0.28, RMSE = 0.054). Using less penalized models - a Lasso Elastic Net (α = 0.5, λ_min = 0.0032) and a pure LASSO (λ_min =0.0041, α = 1) - additional metrics (intra, diff, extraMD, MK) were retained, increasing the explained variance to 61–64% (R²ₗₒₒcᵥ = 0.61–0.64, RMSE = 0.0295–0.032, respectively).\nB)\nWe thank the reviewer for highlighting this point. We acknowledge that partial volume effects can influence diffusion MRI measures, particularly in cortical gray matter. In our study, we minimized these effects by restricting analyses to cortical regions defined by our RS parcellation atlas and averaging metrics across all voxels within each region, reducing the influence of outlier voxels. While residual partial-volume effects cannot be fully eliminated, this approach is consistent with current best practices in neonatal diffusion MRI studies. Importantly, the pipeline is fully automated, and any systematic partial volume effects would be consistent 10 across subjects, allowing valid longitudinal comparisons and between-group analyses. We have also highlighted this limitation in the manuscript (section 4).\n3. The relatively small sample size, combined with the lack of validation on publicly available datasets, raises concerns about the generalizability of the findings. The results would be more convincing with rigorous validation. Although the authors note the difficulty of acquiring longitudinal preterm participants, the Developing Human Connectome Project (dHCP) contains a large cohort of both preterm and term infants, including longitudinal scans. Leveraging this public dataset would greatly strengthen the conclusions.\nWe thank the reviewer for raising the point regarding validation on publicly available datasets. We carefully considered the Developing Human Connectome Project (dHCP) data; however, according to the metadata from dHCP-2 and dHCP-3 releases, only 18 preterm subjects have longitudinal scans within the specific 324/7 – 343/7 wGA window evaluated in our study. From these, only 9 have a second scan within the TEA evaluated window 385/7 – 411/7. From the 9 subjects with data acquired within these developmental windows, only 8 subjects had both diffusion and RS-fMRI data acquired at both times points. Furthermore, these numbers refer to acquired data before preprocessing. Very likely, not all the available subjects will have both raw diffusion and resting-state fMRI data that would pass quality control and be analyzable.\nGiven this extremely small number, validation using the dHCP cohort would be statistically underpowered and could not provide a meaningful assessment of generalizability. Importantly, our cohort is highly homogeneous, comprising healthy very preterm infants (born before 32 wGA) without major neonatal complications, scanned within a tightly defined developmental window, which increases the reliability of our findings and reduces potential confounding effects.\nTherefore, we have chosen not to follow the Reviewer’s suggestion to leverage the available dHCP data. However, we acknowledge that the inclusion of a validation set would strengthen our conclusion and have included a statement to this effect, as highlighted, in Section “Limitations”: Future studies with a larger, longitudinal, age-matched cohort will be important to strengthen the external validity of our findings.\nTable R2. Summary of the dHCP preterm longitudinal data, grouped per scan age and sex.\ndHCP data - Preterm infants N_participants N_female N_male\nPreterm with longitudinal data 84 41 43\nWith First scan (324/7-343/7wGA) 18 10 8\nWith Second scan (385/7- 35 17 18\n411/7wGA)\nWith Both ranges met 9 6 3\nRanges met + both Diffusion 8 6 2\nand fMRI data acquired\n4. While the authors chose 9 microstructure measures from diffusion MRI, BOLD variability was used as the only measure derived from resting state fMRI. Other relatively more widely adopted measures from resting-state fMRI such as amplitude of low frequency fluctuations (ALFF) or regional homogeneity (ReHo) should be included to test the findings or the authors should justify exclusion of these other more widely used measures from resting state fMRI. Lacking these measures from resting state fMRI weakens the functional analysis and limits the ability to fully address how microstructural and biological development support functional maturation.\nWe thank the reviewer for their comment. As noted, BOLD variability has been operationalized using slightly different metrics, including standard deviation (SD), as well as the closely related amplitude of low-frequency fluctuations (ALFF/fALFF). Both these measures aim to capture the amplitude of spontaneous moment-to-moment regional BOLD signal oscillations. In contrast, ReHo (Regional homogeneity) measures local functional connectivity - specifically, the synchronous activity of a voxel with its nearest neighbors - and is not a measure of BOLD variability. This local functional synchronous activity (local correlated-ness) is measured through Kendall correlation between RS-fMRI BOLD signal at the individual voxel (volume unit) level and its nearest neighbors (local connectivity), which is very different from BOLD signal variability.\nOur study aimed to evaluate specifically BOLD variability and its underpinnings. More specifically, how the observed changes in BOLD variability align with regional measures of cortical microstructure, and the patterns of gene expressions that could explain the biological correlates underlying the observed microstructural 12 and functional changes. For this reason, calculating the ReHo is not within the scope of this work.\nHowever, to address the reviewer’s suggestion, we have additionally performed the same longitudinal and cross-sectional analysis using ALFF and fALFF.\nBoth ALFF/fALFF and BOLD SD capture similar aspects of BOLD signal fluctuation over time (Lim et al., 2021; Millar et al., 2020). While ALFF focuses on a specific low-frequency band (0.01–0.08 Hz), SD summarizes overall temporal variance. fALFF further normalizes the amplitude of spontaneous BOLD fluctuations within the low-frequency band (0.01–0.08 Hz) by the total spectral power across all frequencies, aiming to reduce the influence of noise and nonspecific global signal contributions. However, some studies suggest that this normalization may reduce reliability relative to ALFF (Jia et al., 2020; Sbaihat et al., 2022). Overall, since the dominant power of the BOLD signal lies within the low-frequency range (0.01–0.1 Hz), all these measures are highly comparable. As shown in Figure S1 and S5, the longitudinal and cross-sectional results obtained using ALFF/fALFF are consistent with those from BOLD SD, supporting the validity of our findings.\nWe chose BOLD variability as our primary measure because it captures moment-to-moment fluctuations across the full frequency spectrum, providing a broad and sensitive index of spontaneous neural dynamics. ALFF and fALFF are restricted to a predefined low-frequency band and are more dependent on preprocessing choices, while BOLD variability offers a more general and assumption-free characterization of spontaneous activity, which may be better suited to capturing the physiological and microstructural changes occurring during late gestation.\nAs requested, for completeness, we performed parallel analyses with ALFF and fALFF, which yielded consistent results. These results have been referenced in the Results and Methods sections, as highlighted in the manuscript, and have been added to the Supplementary Material.\nResults, Section “Longitudinal cortical BOLD variability and diffusion microstructural changes during early preterm brain development”:\nComplementary analysis using ALFF and fALFF to estimate BOLD variability yielded consistent longitudinal patterns across the same regions (see Supplementary Figure S1).\nResults, Section “Effect of preterm birth on BOLD variability and cortical microstructure at TEA”:\nComplementary analysis using ALFF and fALFF to estimate BOLD variability yielded consistent cross-sectional patterns across the same regions (see Supplementary Figure S5).\nMethods, Section “Regional cortical BOLD variability estimation”:\nIn addition, ALFF (Amplitude of Low-Frequency Fluctuations) was computed by transforming each regional time series into the frequency domain using the Fast Fourier Transform to obtain the power spectrum (Feng et al., 2024; Huang et al., 2020; Yan et al., 2021). The square root of the power spectrum within the 0.01–0.08 Hz (low) frequency band was averaged to obtain regional ALFF values. 13 Regional fALFF (fractional ALFF) was then calculated as the ratio of the low-frequency (0.01–0.08 Hz) ALFF to the total ALFF across all frequencies of the power spectrum. These measures were used in parallel to confirm that the results of BOLD SD analyses of spontaneous BOLD fluctuations were consistent across other different metrics estimating BOLD variability.\nWe selected BOLD variability (BOLD SD) as the primary measure because it captures moment-to-moment fluctuations across the full frequency spectrum, providing a broad and sensitive index of spontaneous neural dynamics. In contrast, ALFF and fALFF are restricted to a predefined low-frequency band and are more dependent on preprocessing choices, whereas BOLD variability provides a more general, assumption-free characterization of spontaneous activity, which is better suited to capturing physiological and microstructural changes during late gestation. ALFF and fALFF results are presented in Supplementary Material.\nFigure S1. Longitudinal BOLD, ALFF and fALFF changes from 33 to 40 weeks. Bar plots illustrate the group-averaged BOLD variability, ALFF and fALFF longitudinal changes in VPT infants from 33 wGA to TEA across networks, with each metric scaled to its own largest absolute mean regional change. Embedded boxplots show the distribution of individual subject changes within each RSN for each metric, with individual values scaled to the largest absolute subject-level change observed across all RSNs for that metric. Significant changes are indicated by asterisks (“*”p<0.05, “**”p<0.01, “***”p<0.001, after FDR correction). PCUN= precuneus, PCC = posterior cingulate cortex, SSM = sensorimotor, VIS = visual, AUD = auditory, THAL = thalamus, PFC = prefrontal cortex.\nMinor Comments:\n5. While the Introduction provides useful background on diffusion models, it does not sufficiently justify why DTI measures (e.g., FA, MD) were applied to study cortical microstructure, given their traditional use in white matter.\nWe thank the reviewer for their comment. Although DTI has traditionally been applied to WM, it has also been applied and successfully used to study GM. As noted in the introduction, DTI does not hold intrinsic diffusivity assumptions, capturing variations in diffusion patterns, reason why it is one of the diffusion models widely accepted to be used both in WM and GM: “Traditional signal representation models of diffusion, such as Diffusion Tensor Imaging (DTI), provide insights into GM organization by measuring the anisotropic diffusion of water molecules within the tissue, without intrinsic diffusivity assumptions, capturing variations in diffusion patterns (Le Bihan et al., 2001)”.\nImportantly, several neonatal and fetal MRI studies have demonstrated that DTI-derived measures are sensitive to early cortical microstructural maturation. In fact, it has been very useful to detect and to state cortical complexification occurring before term age, during the third trimester of pregnancy. During this period there 15 is an important increase in dendritic arborization, and disappearance of the radial glia, leading to a more complex multidirectional and less radial diffusion environment (Bystron et al., 2008; Cameron & Rakic, 1991; He et al., 2020; Marin-Padilla, 1992; Rakic, 2003), which leads to the well described decrease in cortical FA, observed when using DTI to study early cortical maturation (Ball et al., 2013; Batalle et al., 2019; Eaton-Rosen et al., 2015; McKinstry et al., 2002). For this reason, we included DTI in our analysis as a standard baseline diffusion model, allowing us to build on a well-established literature on cortical DTI findings in the developing brain, and we complemented it with the more advanced SMT and DKI models.\n6. Since the authors used three different diffusion models, it would be beneficial to specify the b-values used in each diffusion model in the Methods.\nThank you for the suggestion. All diffusion models (SMT, DKI and DTI) were estimated using the same multi-shell acquisition. To be more explicit, we have added this to Methods, Section “Diffusion cortical microstructural maps”:\nAfter pre-processing, MSDI data was analyzed using three different methodological approaches, each applied to the same set of b-values (b=200, b=1000, b=2000 s/mm2).\nIn the same section, we also state: “For both DTI and DKI metrics, we used the code publicly available from DIPY library (https://docs.dipy.org/stable/interfaces/reconstruction_flow.html), applying the DKI model to the full MSDI data.” This approach avoids inconsistencies that may arise when DTI is fit only to a low-b subset and improves tensor stability by incorporating the additional high-b information.\n7. The definition of “clusters” is not clear. In BOLD variability, “clusters” appear to denote groups of regions showing significant increases between two timepoints. But in diffusion MRI, k-means clustering (k=3) was used without justification of why three clusters were chosen. To avoid confusion, the authors could consider using “groups” for the BOLD variability groups of regions and provide validation (e.g., elbow method, silhouette score) for the k-means choice of clusters.\nWe thank the reviewer for their suggestion. We have changed the definition of “clusters” to “groups” when referring to the regions with or without significant BOLD variability changes.\nIn addition to this grouping based on BOLD variability changes, we wanted to assess if the longitudinal microstructural maturational changes would result in a similar regional pattern grouping (independently of BOLD variability changes). For that reason, we have conducted an additional clustering analysis, as described in Methods, Section “MRI data analysis” and Supplementary Material Figure S3.: “Microstructural diffusivities were clustered using the R-package and function ConsensusClusterPlus (k-means base algorithm and Euclidean distance) (Wilkerson & Hayes, 2010), to identify stable clusters of regional diffusivities. The optimal number of clusters was determined using the elbow method, based on the within-cluster sum of squares (WCSS).” We selected k = 3 as the optimal number 16 of clusters using the elbow method, which showed a clear inflection point at this value. Beyond k = 3, the reduction in within-cluster sum of squares (WCSS) diminished substantially, indicating limited gain in compactness with additional clusters. We have added Figure S2 a), which displays the elbow method.\nWe acknowledge that the consensus clustering analysis results have changed slightly following revision. This is because we now used the group-level region × metric matrix (representing the relative spatial distribution of longitudinal changes across cortical resting-state network regions) preserving the directionality of effects, as explained in Methods section. For clustering, retaining the sign of metric changes is important, as positive versus negative changes reflect distinct patterns of microstructural maturation across regions. This contrasts with Lasso-based regression, which focuses on the magnitude of change, and for which the directionality of effects is not informative for feature selection. Importantly, these updates do not alter the main findings or overall interpretation.\nAs explained in Figure S2 paragraph: “Three distinct regional microstructural maturational patterns were identified, aligning with the central-to-peripheral and posterior-to-anterior known gradients in brain maturation, in agreement with the brain myelination order of Kinney (Kinney et al., 1988). Microstructural cluster 2 aligns with BOLD SD group 1, while Microstructural cluster 3 is more similar to BOLD SD group 2. The thalamus emerged as a separate cluster, reflecting its distinct microstructural profile.” This microstructural spatial distribution, although not perfectly identical, is similar to the regional patterns identified by BOLD variability changes, highlighting that the primary sensory-motor and proto-DMN regions exhibit both functional and microstructural changes distinct from those of other brain regions”.\nFigure S2. Results from the consensus clustering of the combined 9 distinct\na) Elbow method showing an inflection point at K = 3, suggesting an optimal clustering solution with three clusters. b) Heatmap depicting the regional maturational patterns from 33- to 40-wGA, with rows and columns representing the RSNs and color intensity in each cell reflecting the composite changes over time of the microstructural metrics. c) Composition of the identified 3 regional clusters of microstructural delta changes (Cluster 1: Thalamus, Cluster 2: PCC, Visual, PCUN, SSM, Auditory, Limbic; Cluster 3: PFC, Paralimbic). d) Brain plots illustrating the 3 distinct regional clusters depicted by the consensus clustering analysis, color-coded according to regional cluster (red = cluster 1, orange = cluster 2, blue = cluster 3).\n17\n8. The authors might want to be cautious about using “trajectory”, as only two timepoints were included in the manuscript.\nWe thank the reviewer for their suggestion and have removed the word “trajectory” from the discussion.\n9. There is a typo in the last sentence of the Introduction: “weather” should be corrected to “whether.”\nWe thank the reviewer for detecting this typo and we have corrected it accordingly.\nReviewer #3 (Remarks to the Author):\nThank you for asking me to review this interesting paper examining neonatal brain development. The work combines MRI data with open access transcriptomic data to produce interesting and plausible biologic insights. The combination of BOLD as a measure of cortical function with transcription data is as far as I know completely novel and of great interest. The combination of MRI with transcription data is an emerging and promising area of science and Ball is the world leader in developing this approach in the developing brain. The paper leverages longitudinal imaging data encoding the 'experiment of nature' that is preterm delivery. Premature birth is a wide-ranging and robust disrupter of brain development, which has been shown (by these authors as well as many others) to affect almost every aspect of brain development and so provides a useful 'positive control' reality check on these results.\n18\n1. The imaging data appropriately collected with reasonably up to date sequences. The cohort is small, and indeed quite a lot smaller than most substantial studies currently being published. The data were collected by the researchers, and they probably trust their own data quality, relying on the longitudinal data to add value to the small cohort. Nevertheless, in the same way that they used open access transcriptomic data, they could have replicated their results using the open access imaging data now freely available, which are an order of magnitude larger than the current dataset, and which some of this group have used in the past.\nWe thank the reviewer for their relevant comment. We carefully chose a longitudinal design for our study because both BOLD-variability and microstructural diffusion metrics are subject-dependent, and single time-point cross-sectional group comparisons are less reliable for capturing developmental changes. Longitudinal within-subject data, collected at well-defined time-points, allows us to more robustly track individual trajectories of cortical maturation, which is particularly critical during this period of rapid developmental changes.\nAs noted in response to Reviewer 2 above, we have considered the Developing Human Connectome Project (dHCP) data; however, according to the metadata from dHCP-2 and dHCP-3 releases, only 18 preterm subjects have longitudinal scans within the specific 324/7 – 343/7 wGA window evaluated in our study. From these, only 9 have a second scan within the TEA evaluated window 385/7 – 411/7. From the 9 subjects with data acquired within these developmental windows, only 8 subjects had both diffusion and RS-fMRI data acquired at both times points. Furthermore, these numbers refer to acquired data before preprocessing. Very likely, not all the available subjects will have both raw diffusion and resting-state fMRI data that would pass quality control and be analyzable. Given this extremely small number, validation using the dHCP cohort would be statistically underpowered and could not provide a meaningful assessment of generalizability. Importantly, our cohort is highly homogeneous, comprising healthy very preterm infants (born before 32 wGA) without major neonatal complications, scanned within a tightly defined developmental window, which increases the reliability of our findings and reduces potential confounding effects.\nTherefore, we have chosen not to leverage the available dHCP data. However, we acknowledge that the inclusion of a validation set would strengthen our conclusion and have included a statement to this effect, as highlighted, in section “Limitations”: Future studies with a larger, longitudinal, age-matched cohort will be important to strengthen the external validity of our findings.\nTable R2. Summary of the dHCP preterm longitudinal data, grouped per scan age and sex.\ndHCP data - Preterm infants N_participants N_female N_male\nPreterm with longitudinal data 84 41 43\nWith Second scan (385/7- 35 17 18\n411/7wGA)\nWith Both ranges met 9 6 3\nRanges met + both Diffusion 8 6 2\nand fMRI data acquired\nTable R3. Subject-level description of the dHCP preterm longitudinal data within the specific GA window evaluated in our study, including sequence acquisition (MSDI and fMRI) availability.\n2. The image analysis is reasonable. It appears that there was no motion correction applied to the scan data, but a reasonable cut-off for exclusion of scans was used to remove corrupted datasets. It is a bit controversial whether this scrubbing of individual data (which disrupts the time course) affects fMRI interpretation, but my view is that it is unlikely to invalidate measures of BOLD variability. Variability is a nice and under-explored measure of cortical function, the results are interesting and as far as I know novel, and the functional locations revealed are interesting. Spherical means is a good method for analysing diffusion data and the results are consistent with previous studies using different approaches.\nWe thank the reviewer for this comment. The measure of BOLD variability used in our study reflects the overall amplitude of signal fluctuations (standard deviation across time) and does not rely on the temporal sequence or frequency characteristics of the signal. Therefore, removing high-motion volumes (“scrubbing”) does not bias or distort the amplitude measure, as it would only minimally affect the overall variance calculation. By contrast, metrics such as the Hurst exponent (HE) or the mean squared successive difference (MSSD) depend on temporal correlations or the temporal ordering of time points. Scrubbing disrupts the time series, which can strongly influence these measures. In this 20 sense, BOLD variability is inherently more robust to motion-related data removal, supporting the reliability of our results despite excluding corrupted scans.\n3. The imaging results are plausible and interesting, with the increases in visual and sensory-motor cortex in late gestation as expected given known developmental pathways, and to is also interesting that they did not find BOLD variability increases in pre-frontal cortex or limbic- it is plausible to interpret this as the frontal cortex not having started its rapid development while the limbic system may have completed its comparable phase. However, it would have been interesting to see more data on temporal cortex which is rapidly developing at this time.\nWe thank the reviewer, and we agree with these important insights. Our resting-state networks (RSNs) atlas was generated from the combination of RSNs identified in our cohort of infants, including both VPT at 33- and 40-wGA, as well as the FT infants, combined into a cortical parcellation atlas. We chose this approach since we preferred to use RSN directly generated from our data during this developmental time, rather than applying an external atlas to our data. Consequently, our analysis is constrained to the RSN that were observed. Notably, we observed an increase in BOLD variability in the auditory network from 33 to 40 weeks' GA. The temporal pole was included within the paralimbic belt network, together with the orbito-frontal cortex, and this network did not show such an increase in BOLD variability. This further supports a different maturational pattern between primary vs paralimbic/associative regions.\n4. The transcriptomic data and analysis are high quality and appropriate. it is helpful to have the three time points in the data even if there are no imaging data to compare to the first period. The gene set analysis is uncontroversial and Lasso is a reasonable approach to this underdetermined dataset. As with the imaging data the results are plausible and it is re-assuring that neuronal ensheathment is most significant relation, which gives confidence in other pathways detected. It is a shame that the analysis does not make use of the remarkable microbrain histologic parcellation that one of the authors has recently published but relies on small patch location for the transcriptomic analysis, but that does not invalidate the results. The gene sets exposed in the analysis are interesting and add to existing knowledge.\nWe thank the reviewer very much for their positive feedback. As discussed above in response to Reviewer 1, despite a growing catalogue of prenatal gene expression studies, including those that span the whole of gestation (e.g.: BrainSpan), or those with multiple cortical regions (e.g. μBrain), none yet exist that offer both spatial and temporal resolution during this time period. As our window of observation was focused on the third trimester, we chose to leverage the temporal resolution of BrainSpan to examine trends in gene expression over the prenatal window. The cortical regions sampled in the BrainSpan data also cover major regions of the neocortex, allowing for a regional analysis over time. The μBrain atlas, while offering more detailed spatial resolution, is based on microarray data from a small window in mid-gestation (15-20weeks), right at the onset of gliogenesis, and, given the rapid structural, cellular and morphological changes during mid- to- late gestation, it is not certain how well this would map onto the imaging observations acquired toward the end of gestation. We anticipate that, with new resources, exploiting both spatial and temporal dimensions will yield interesting future research directions in this area, including the structural and functional maturation of the cortex during gestation."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I thank the authors for carefully addressing my comments. I believe the updated version of this manuscript is much improved and will be of great interest to the field.", "reviewer_2": "It is well appreciated that the authors have done a great job in revising this manuscript. The authors have addressed previous comments by incorporating additional and relevant references on BOLD variability, conducting further analyses using other BOLD signal measures, providing more detailed methodological descriptions, and discussing the limitations related to the relatively limited dataset. The manuscript has been improved significantly after revision. Just a couple of minor comments remain as follows.\n1. It would be helpful to clarify in the Methods section whether all features included in the LASSO regression were standardized prior to model fitting. Given the different scales of the diffusion-derived measures, explicit reporting of normalization is important. In addition, several of the microstructural metrics are highly correlated; therefore, the selection of MD as a predictor should be interpreted with caution. The authors may consider toning down language suggesting MD as the “best” predictor and instead emphasize that LASSO identifies representative features from correlated sets.\n2. The Discussion cites prior studies applying NODDI to cortical gray matter. Since NODDI was established for white matter modeling and does not account for somas, making its application to cortical gray matter problematic. Previous studies tailored for cortical gray matter microstructure (e.g. PMID: 32289460; PMID: 41486548) have repeatedly found that NODDI might not be used appropriately for quantifying cortical microstructure. Considering these specific works, the authors should briefly acknowledge the methodological limitations and interpret NODDI-related findings with caution.", "reviewer_3": "The authors have addressed my questions."}, "author_response": "We thank sincerely the reviewers for their positive feedback. Detailed responses are provided below in blue.\nReviewer #1 (Remarks to the Author):\nI thank the authors for carefully addressing my comments. I believe the updated version of this manuscript is much improved and will be of great interest to the field.\nReviewer #2 (Remarks to the Author):\nIt is well appreciated that the authors have done a great job in revising this manuscript. The authors have addressed previous comments by incorporating additional and relevant references on BOLD variability, conducting further analyses using other BOLD signal measures, providing more detailed methodological descriptions, and discussing the limitations related to the relatively limited dataset. The manuscript has been improved significantly after revision. Just a couple of minor comments remain as follows.\n1. It would be helpful to clarify in the Methods section whether all features included in the LASSO regression were standardized prior to model fitting. Given the different scales of the diffusion-derived measures, explicit reporting of normalization is important. In addition, several of the microstructural metrics are highly correlated; therefore, the selection of MD as a predictor should be interpreted with caution. The authors may consider toning down language suggesting MD as the “best” predictor and instead emphasize that LASSO identifies representative features from correlated sets.\nWe thank the reviewer for their suggested, we have now clarified in the methods section that the features included in Lasso regression were normalized prior to model fitting. “For Lasso-based regression modelling, the same normalized group-level region × metric matrix was used (mean regional changes were first computed for each cortical resting-state network region and then normalized per metric by dividing by the largest absolute mean regional change), but expressed as absolute values, in order to focus on the magnitude of change, irrespective of the direction (…).”\nWe also had already mention it in the following section: “The analysis was conducted on the group average 9 region x 10 metric matrix (n = 9 regions), where BOLD variability change served as the dependent variable and nine microstructural metrics (extraMD, extrTrans, diff, intra, MD, FA, MK, AK, RK) served as predictors. All metrics were standardized prior to modelling, as previously described.”\nWe also agree with the reviewer suggestion regarding interpretability of Lasso results and we have removed the sentence: “(…) with mean diffusivity (MD) showing the most robust association with longitudinal changes in BOLD variability” from the discussion section. Now it can be read as follows. “Diffusion- and kurtosis-derived microstructural metrics emerged as strong contributors. Under strict penalization, a limited set of predictors (MD, RK, AK and FA) was retained, explaining up to 28% of the observed changes in BOLD variability. Relaxing model sparsity increased the number of contributing microstructural metrics, including intra, diff, extraMD, MK, and substantially improved predictive performance up to 61–64%, suggesting that BOLD variability reflects a distributed and correlated microstructural maturation process. Together, these findings support an important role for cortical microstructural maturation in shaping functional variability.”\n2. The Discussion cites prior studies applying NODDI to cortical gray matter. Since NODDI was established for white matter modeling and does not account for somas, making its application to cortical gray matter problematic. Previous studies tailored for cortical gray matter microstructure (e.g. PMID: 32289460; PMID: 41486548) have repeatedly found that NODDI might not be used appropriately for quantifying cortical microstructure. Considering these specific works, the authors should briefly acknowledge the methodological limitations and interpret NODDI-related findings with caution.\nWe thank the reviewer for their suggestion. We had already acknowledged the limitations of NODDI in cortical gray matter in the Introduction, when justifying our choice of SMT and DKI models. The paper cited in the Discussion primarily uses other microstructural models and reports that “primary sensory and motor cortices exhibited microstructural features indicating more advanced maturation, such as higher T1/T2 contrast, higher intracellular volume fraction, and lower mean diffusivity, compared with higher-order frontal and parietal regions.” These microstructural findings are broader than the individual NODDI intra-cellular volume fraction and retaining the reference in this context does not invalidate their results, despite the known limitations of NODDI. For this reason, we have kept this discussion point while already noting NODDI’s limitations in the Introduction.\nReviewer #3 (Remarks to the Author):\nThe authors have addressed my questions."}]

Decision letter