We acquired and analyzed two sets of research MRI data, one from full term healthy neonates and one from healthy 8-year-old children. The demographic information of the research subjects is summarized in Table 1, and there was no difference in sex between the two groups (p > 0.05). Informed consent was obtained from all subjects following Institutional Review Board approved study protocols at Arkansas Children’s Hospital. Anonymized data were used for analysis.

The neonates and 8-year-old children’s MRI data were both acquired using a 3.0 T Prisma (Siemens healthcare) scanner. The neonatal scans were performed when they were at natural sleep with no sedation. The 8-year-old children were instructed to close their eyes during the RS-fMRI data acquisition. The MRI protocol included a MPRAGE 3D T1-weighted scan with TR 2400ms, TE 2.24ms, 1 average, voxel size 0.8*0.8*0.8mm³, 8° flip angle, turbo factor 256, and sagittal slices covering the entire brain, and a RS-fMRI pulse sequence with multi-band EPI with TR 800ms, TE 37ms, phase-encoding AP direction then repeated in PA direction, multi-band factor 8, matrix size 104*104, voxel size 2*2*2 mm³, and 72 slices. In addition, a short field map sequence using the same geometry was performed for both AP/PA phase encoding, to correct for field inhomogeneity associated distortion in imaging post-processing.

The 8-year-old children’s brain structural and functional images were preprocessed by the fMRIPrep toolbox (Esteban et al., 2019), which includes tools from FSL (Jenkinson et al., 2012), AFNI (Cox, 1996), and ANTs (Avants et al., 2009). The pipeline in our study includes (1) slice timing correction; (2) correction for susceptibility distortions induced by magnetic field inhomogeneity and multi-band EPI; (3) realignment of all volumes to a selected reference volume; (4) co-registration of the functional data to the structural image; (5) normalization to the MNI standard space, data resampling to 2*2*2 mm³ isotropic voxels; (6) functional images’ segmentation into gray matter and white matter; and (7) smoothing the segmented white matter with a 4 mm FWHM Gaussian kernel.

Since the fMRIPrep toolbox by default lacks an age-appropriate template for the neonatal brain, we applied in-house code base on FSL, AFNI and ANTs, same as those used in fMRIPrep, to perform data preprocessing on neonates (Jenkinson et al., 2012; Cox, 1996; Avants et al., 2009). It should be noted that all steps were the same as those used in fMRIPrep and the only difference between processing of the neonatal and 8-year-old data was that the neonatal data were registered to the UNC neonatal template (Shi et al., 2011) whereas data for 8-year-old children were registered to the MNI space. And the white and gray matter masks from UNC were used to segment the white and gray matter data for the neonates.

Head motion of all subjects was checked, and only subjects with RS-fMRI data with minimal to slight head motion (<3.0 mm in translation and < 3° in rotation for all imaging volumes for 8-year-old children; <2.0 mm in translation and < 2° in rotation for all imaging volumes for neonates) were kept in the study. After the head motion evaluation, 3 subjects from the 8-year-old group were excluded and 38 subjects were used for further analysis; and 3 subjects from the neonatal group were excluded and 69 subjects were used for further analysis.

We first decomposed the two preprocessed BOLD data separately using standard group-level ICA to identify statistically ICs, as implemented in the GIFT toolbox.¹ The data were decomposed using the Infomax ICA algorithm, and the ICASSO algorithm was repeated 20 times to increase the stability of ICs. Subject-specific spatial maps and time courses were estimated using the GICA back-reconstruction method, similar to that performed in previous research on whole-brain connectivity in the resting state (Allen et al., 2014). FMRI data of white matter were decomposed into 16 and 20 independent components for the neonates and 8-year-old children, respectively. Eight year old children have a more mature brain network development compared to neonates, and their BOLD signal were therefore decomposed to a larger number of independent components, consistent with previous studies (Gao et al., 2015b). A similar processing scheme was also applied to the gray matter fMRI data for references.

The group-average spatial maps obtained from the GIFT output were visually inspected, with each component manually labeled as either a functional network or noise. Following the method of functional network identification in gray matter ICA decomposition by Du et al. (2023), we organized the spatial maps from the white matter ICA decomposition into specific functional networks. To determine whether a component constituted a meaningful functional network and should be included in the analysis, we removed components located on the scalp or blood vessels, compared each component to previous research findings (Gao et al., 2017; Gilmore et al., 2018; Allen et al., 2014; Allen et al., 2011; Seeley, 2019; Du et al., 2020), and identified white matter functional networks according to their anatomical roles and the JHU atlas (Du et al., 2023; Allen et al., 2011).

For the comparison of functional networks between the two age groups, since the functional brain networks of neonates are still developing and not fully mature (Huang Z. et al., 2020; Gao et al., 2015b; Gao et al., 2015a), it may be challenging to establish a one-to-one correspondence between the independent components obtained from neonates and those obtained for 8-year-old children. We focused on optic radiation (OR), cerebral peduncle (CP) and anterior corona radiation (ACR) functional networks in both datasets, which are consistently delineated in both groups (Figure 1) and are essential for the visual, motor, and cognitive functions. For gray matter, we identified the primary functional networks including the visual network (VN), the sensorimotor network (SMN), and the auditory network (AN), as well as high-order networks including the default mode network (DMN), the salience network (SAN), and the frontoparietal networks (FPN) in both datasets (Supplementary Figure S1).

To remove the influence of irrelevant signals, we conducted the following post-processing on all time courses (TC): transforming TC to a Fisher’s Z-score, regressing six head motion parameters, and removing linear trends. This analytical framework is systematically applied to both white matter and gray matter fMRI data sets respectively, ensuring a comprehensive understanding of FC across different brain tissue types.

Consequently, we analyzed the differences in intra-network and inter-network FCs among the 3 white matter functional networks and 6 gray matter functional networks between the two age groups. As an additional measure, we also examined the differences in both intra-network FC and inter-network FC for all meaningful functional networks (as presented in Figure 1) combined, regardless of whether they were consistently detected in both age groups.

For each functional network, we assessed intra-network FC by computing the mean Pearson correlation coefficient across the voxel TC within the network, excluding diagonal elements. For each participant, we also calculated the average intra-network FC for all functional networks to derive the intra-network FC. The main objective of intra-network FC is to assess the synchronous activity within the functional network. Positive correlations reflect the temporal consistency of activities among different voxels within the network, which is crucial for measuring the network’s homogeneity. Negative correlations, on the other hand, might indicate functional differentiation within the network, which is not aligned with our purpose of assessing homogeneity. In brain functional studies, positive correlations are often viewed as indicators of functional connectivity, representing regions that are either simultaneously activated or inhibited. Negative correlations are less commonly used to denote such synchronous activity, especially when evaluating connectivity within a single network. Therefore, for intra-network FC, we only considered positive correlations within the network.

For each participant, we computed the inter-network FC by evaluating the mean Pearson correlation coefficient matrix among all pairs of the identified functional networks, again omitting diagonal elements. This process provides us with distinct inter-network FC values for each individual. For inter-network FC, we include both positive and negative correlations between networks. The goal of inter-network FC is to assess the interactions between different functional networks. Positive correlations indicate cooperative activity between networks, while negative correlations represent a complementary or inhibitory relationship. Both types of correlations are biologically significant in the context of inter-network interactions. Considering the complex interaction patterns between different networks, including both positive and negative correlations ensures a more comprehensive evaluation of inter-network functional connectivity. Omitting negative correlations would overlook crucial information about the nature of these interactions. Averaging both positive and negative correlations to evaluate overall inter-network functional connectivity is commonly used in clinical and cognitive neuroscience research, especially in studies that require a comprehensive assessment of network connectivity strength or changes in connectivity patterns (Berman et al., 2016; Sahoo et al., 2015; Song et al., 2015). Therefore, we considered both positive and negative inter-network FC.

The amplitude of low-frequency fluctuation (ALFF) of the resting-state fMRI signal has been suggested to reflect the intensity of regional spontaneous brain activity (Zang et al., 2007). Due to its sensitivity to noise, the fALFF was subsequently proposed to study brain function (Zou et al., 2008). Here we calculated the ratio of power spectrum at 11 frequency bands to that of the frequency range (0.01 ~ 0.15 Hz), which denote the power of BOLD signal fluctuation. These 11 frequency bands starting from 0.01 Hz to 0.15 Hz, each frequency band is 0.04 Hz in length, with a step size of 0.01 Hz. For each independent component of each subject, TC was extracted to compute fALFF at the 11 frequency bands. Similar to the process used in calculating FC, we also applied the following post-processing to TC: Fisher-Z score transformation, head-motion regression, and detrending before computing fALFF.

The identified white matter functional networks for neonates are as follows: OR (Optic Radiation), GCC (Genu of Corpus Callosum), CP (Cerebral Peduncle), SCR (Superior Corona Radiata), ACR (Anterior Corona Radiata), PCR (Posterior Corona Radiata), and SS (Sagittal Stratum), as shown in Figure 1a. The identified white matter functional networks for 8-year-old children are as follows: OR (Optic Radiation), SCC (Splenium of Corpus Callosum), CP (Cerebral Peduncle), SCP (Superior Cerebellar Peduncle), ACR (Anterior Corona Radiata), MCP (Middle Cerebellar Peduncle), PLIC (Posterior Limb of Internal Capsule), RPIC (Retrolenticular Part of Internal Capsule), and UF (Uncinate Fasciculus), as shown in Figure 1b. Possibly due to the immaturity of functional networks in neonates and the much lower BOLD signal in white matter, only three networks: OR, CP and ACR were consistently identified (the first three rows of the left column in Figures 1a,b) and were used for group comparison.

To investigate how functional networks evolve from neonates to 8-year-old children, we compared the intra-network FC of the above three networks, with the corresponding results presented in Figure 2, which shows that the intra-network FC within OR, CP, and ACR were lower (all p values <0.05, FDR corrected) when compare 8-year-old children to neonates.

To explore the differences in intra-network FC of all white brain networks between neonates and 8-year-old children from a holistic perspective, we further examined the overall intra-network FC within white matter for all functional networks derived from the ICA algorithm (Figure 1) for the two age groups, and the results are shown in Figure 3. It can be also seen that the intra-network FC was lower in 8-year-old children compared to that in neonates. These results suggest that with an increase in age, the intra-network FC within the same functional networks decreases.

Next, we compared the inter-network FC of the three networks between neonates and 8-year-olds, with the corresponding results presented in Figure 4, which indicated that the inter-network FC between CP and ACR were higher when compare 8-year-old children to neonates (p < 0.05), while the inter-network FC between OR and ACR were similar between the two age groups (p = 0.08). Here we mainly measure the inter-networks between basic and advanced networks, so we did not present the inter-networks FC between OR and CP (which are both basic brain functional networks).

Similarly, we also examined the overall inter-network FC within white matter for all functional networks using the ICA algorithm (Figure 1) and compared them between the two age groups, and the results are shown in Figure 5. It can be also seen that the inter-network FC was higher in 8-year-old children compared to those in neonates. These results suggest that with an increase in age, the inter-network FC between the different functional networks increases.

Firstly, we calculated the fALFF of each component for each participant (excluding noise components). Then, we averaged the fALFF of all components for each participant to obtain the fALFF of that participant in a specific frequency band. Finally, we compared the fALFF across different frequency bands between neonates and 8-year-old children. Multiple comparison corrections (FDR) were applied to account for the 11 frequency bands. The comparison of fALFF in white matter functional networks is presented in Figure 6. Specifically, in the low-frequency range (0.01–0.05 Hz, 0.02–0.06 Hz, 0.03–0.07 Hz, 0.04–0.08 Hz), the fALFF for the 8-year-old children were lower than those for the neonates (p < 0.05, FDR corrected); and in higher frequency bands (0.07–0.11 Hz, 0.08–0.12 Hz, 0.09–0.13 Hz, 0.10–0.14 Hz, 0.11–0.15 Hz), the 8-year-old children had higher fALFF than the neonates (p < 0.05, FDR corrected).

We also found similar results in the gray matter functional network: the intra-network FC of 8-year-old children is lower than that of newborns, while the inter-network FC is higher than that of newborns. Additionally, with increasing age, the distribution of fALFF shifts toward higher frequencies. Detailed results can be found in the Supplementary material.