Resting state functional networks in 1-to-3-year-old typically developing children.

Brain functional networks undergo substantial development and refinement during the first years of life. Yet, the maturational pathways of functional network development remain poorly understood. Using resting-state fMRI data acquired during natural sleep from 24 typically developing toddlers, ages 1.5-3.5 years, we aimed to examine the large-scale resting-state functional networks and their relationship with age and developmental skills. Specifically, two network organization indices reflecting network connectivity and spatial variability were derived. Our results revealed that reduced spatial variability or increased network homogeneity in one of the default mode network components was associated with age, with older children displaying less spatially variable posterior DMN subcomponent, consistent with the notion of increased spatial and functional specialization. Further, greater network homogeneity in higher-order functional networks, including the posterior default mode, salience, and language networks, was associated with more advanced developmental skills measured with a standardized assessment of early learning, regardless of age. These results not only improve our understanding of brain functional network development during toddler years, but also inform the relationship between brain network organization and emerging cognitive and behavioral skills.

Introduction

A vital goal of developmental neuroscience is to characterize brain development during infancy and early childhood, a key period of postnatal brain maturation (Tau and Peterson, 2010), peak neuroplasticity (Kolb and Gibb, 2011), and a critical window for emergence of cognitive and behavioral skills (Johnson et al., 2008). Atypical brain maturation has been implicated in various neurodevelopmental and neuropsychiatric disorders, such as autism spectrum disorder and attention deficit hyperactivity disorder (Bernhardt et al., 2017; Posner et al., 2014), major depressive disorder (Mulders et al., 2015), and schizophrenia (van den Heuvel and Fornito, 2014). Yet, our current understanding of brain development early in life is still limited, in part due to the challenges of acquiring in-vivo brain imaging data in very young children.

A number of cross-sectional and longitudinal studies have begun to characterize the maturational pathways of brain structure and functional networks in the first years of life. Generally, these studies have indicated that brain structure and functional networks are already in place by the second year of life, undergoing reorganization and fine-tuning beyond two years of age. Cortical and subcortical gray matter volumes, in particular, grow most rapidly in the first year of life and peak around age 2 years (Gilmore et al., 2012; Haynes et al., 2020; Knickmeyer et al., 2008; Pfefferbaum et al., 1994). This rapid gray matter volume growth likely reflects continuous synapse formation in the first years of life, albeit with variable rates across different brain regions (Huttenlocher and Dabholkar, 1997). Morphometrically, beyond the volumetric growth, maturational changes in brain structure are indexed with dramatic increases in cortical thickness and cortical folding (gyrification) and expansion of cortical surface area observed in the first years of life (Brown and Jernigan, 2012; Li et al., 2013, 2014; Lyall et al., 2015). The rate of these developmental processes varies across the brain, with certain regions growing more rapidly than others (Tau and Peterson, 2010). Unlike this macro-scale cortical structural maturation which peaks around 2 years of age, brain functional networks continue to develop well beyond infancy and into later childhood, especially those implicated in higher-order cognitive functions (Gao et al., 2015a; Hoff et al., 2013).

Current understanding of functional network development has been primarily informed by resting-state functional connectivity, which measures temporal correlations (or statistical dependence) of spontaneous, low-frequency fluctuations of the blood oxygen level-dependent (BOLD) signal in different brain regions, using functional magnetic resonance imaging (fMRI). A common approach for identifying resting-state functional networks (RFNs) is independent component analysis (ICA), a data-driven approach that statistically decomposes fMRI data into a set of spatial components with maximally independent time courses (Beckmann, 2012), with many components showing close correspondence to known neuroanatomy and canonical brain functional activation maps (Smith et al., 2009). Based on their spatial and temporal characteristics, independent components are classified as RFNs (e.g., visual network spanning regions in the occipital visual cortex, motor network spanning the motor strip, etc.) or non-neuronal noise or artifact components (e.g., generated by motion during scanning, or corresponding to physiological variables, such as breathing or pulsation in large blood vessels) (Griffanti et al., 2017).

Utilizing ICA approach, distinct RFNs have been reliably described in adults, including primary, unimodal sensory networks (i.e., visual, auditory, and sensorimotor networks) and higher-order or supra-modal networks, such as self-reference (default mode), language, attention, and executive networks (Damoiseaux et al., 2006; Smith et al., 2009). Primary sensory networks typically exhibit more localized connections, often limited to homologous regions in each hemisphere. In contrast, higher-order networks often encompass spatially distal nodes distributed across multiple brain regions and lobes.

Comparable functional brain networks have also been reliably identified early in life in young children (Fransson et al., 2007; Gao et al., 2015a,b). However, while primary sensory networks undergo subtle refinement and strengthening over the first two years of life, and substantially resemble adult topology by age two (Gao et al., 2015a; Lin et al., 2008), higher-order, supra-modal functional networks are far from the adult-like organization in the first postnatal years, and undergo prolonged maturation over the first decades of life (Dosenbach et al., 2010; Fair et al., 2008; Gao et al., 2015a). A longitudinal study of infants scanned in the first postnatal year at three-month intervals has demonstrated this maturational sequence of functional network connectivity across the first year of life (Gao et al., 2015a). The findings revealed that primary sensorimotor, auditory, and visual networks resemble adult-like topologies soon after birth and exhibit minimal topological changes during the first year, suggesting that they are the earliest developing functional networks. In contrast, higher-order networks such as attention, default mode, and frontoparietal executive networks scarcely resemble their adult topologies at birth and undergo age-dependent increases in connectivity, with some networks such as default mode (but not others) reaching adult-like topology by the end of the first year. Critically, this maturational sequence of RFNs parallels the order in which behavioral and cognitive skills emerge and develop in early childhood (Johnson, 2001).

Since the seminal studies by Gao et al. describing the development of RFNs in infants, the maturational trajectories of RFNs in toddlerhood and preschool age (1–4 years) remain understudied (likely as a result of unique practical and methodological challenges associated with obtaining high quality MRI data in young children, necessitating adapted data acquisition procedures and consideration of age-specific and developmentally-appropriate data processing and analytic approaches; Turesky et al., 2021). Only a few studies have examined brain activation or functional connectivity in toddlers and preschoolers, with a particular focus on speech perception and language circuitry, either through passive auditory stimulation or resting-state during sleep (Hutton et al., 2015; Redcay et al., 2008; Xiao et al., 2016). One study examined local and global functional connectivity features (i.e., amplitude of low frequency fluctuation, regional homogeneity, and eigenvector centrality mapping) in awake 2.5–6 years old children during passive viewing fMRI (Long et al., 2017). The results revealed age-related increase in local and global connectivity in regions within the default mode and frontoparietal networks, and age-related shift from more global to local connectivity in the superior parietal and fusiform gyri, and a local-to-global shift in the superior temporal area, suggesting that functional connectivity undergoes substantial reorganization during toddler and preschool years. Finally, a recent study by Bruchhage et al. (2020) also reported age-related effects in functional network connectivity and links with developmental skills in a cross-sectional cohort spanning a relatively broad age range from infancy to early school age (3 months to 6 years). Using predefined regions of interest derived from adult datasets, this study identified a cross-age shift in functional connectivity toward networks involved in higher-order cognitive processes (e.g., default mode, attention, and salience networks). The results also revealed distinct functional connectivity patterns associated with visual, motor, and language skills measured with behavioral assessments, with developing motor and language skills – independent of age – linked with overlapping connectivity patterns in the default mode, visual, salience, and dorsal attention networks. Outside of these singular studies, however, the evidence on development and organization of large-scale RFNs – and the links between functional network connectivity and emerging cognitive skills – during the toddler years remains scarce.

In the current study, we aimed at examining the brain functional networks in a sample of 1-to-3-year-old typically developing children, using resting-state fMRI acquired during natural nocturnal sleep. We investigated age-related effects for two indices of functional network organization, within-network connectivity and within-network spatial variability, in ICA-derived RFNs. We predicted that, given the maturational sequence of functional networks previously reported in infancy (Gao et al., 2015a), higher-order functional networks would exhibit positive associations with age in within-network connectivity and negative associations with age in within-network spatial variability, with no such associations expected in the primary unimodal network (e.g., visual, sensorimotor networks). In addition, we predicted that network organization indices in higher-order networks would be associated with more favorable developmental outcomes.

Methods and materials

Participants

This study includes data from 24 typically developing toddlers and preschoolers between 1.5 and 3.5 years of age (between 539 and 1249 days old, corrected to a 40-week gestation; see Table 1 for demographic information) enrolled in the San Diego State University (SDSU) Toddler MRI Study. Children were recruited from the community, including Early Head Start and other early childhood education programs, and print and social media advertisements, in and around San Diego County. Participants were excluded for prematurity (<36 weeks of gestation), family history (in first degree relatives) of heritable neuropsychiatric or neurological disorders, personal or family history of neurodevelopmental disorders, such as autism spectrum disorder and intellectual disability, and history of perinatal CNS infection or gross CNS injury, non-febrile seizures, and contraindications for MRI. The research protocol was approved by institutional review boards of SDSU and University of California San Diego (UCSD). Written informed consent was obtained from the caregivers.

Table 1

Participant Characteristics and Demographic Information.

	N = 24
Gender (M/F)	14/10
Ethnicity (Hispanic/Non-Hispanic)	5/19
Race (White/Asian/More-than-one/Black/Unknown)	16/4/2/1/1
Delivery method (Vaginal/C-section)	16/8
Exposure to more than one language at home (Yes/No)	8/16
Maternal education level (%)
High school or some college credit, but < 1 year	8%
Associate degree	4%
Bachelor’s degree	42%
Master’s degree	38%
Professional degree (MD, PhD, JD)	8%

Note: M = male; F = female; RMSD = root mean squared displacement.

CBCL data are missing for one subject.

Mean RMSD and # of volumes censored were calculated across two fMRI runs.

Developmental assessment

Participants’ developmental skills were assessed with the Mullen Scales of Early Learning (MSEL; Mullen, 1995), a clinician-administered standardized assessment of cognitive, language, and motor development. The Vineland Adaptive Behavior Scales, Second Edition, Survey Form (VABS-2; Sparrow et al., 2005), a semi-structured caregiver interview, was administered to parents to assess the child’s adaptive functioning skills and aid with ruling out any developmental challenges. Parents also completed the Social Communication Questionnaire (SCQ, Current form) (Lord and Rutter, 2003), a screener for autism spectrum disorder, with no participants exceeding the cutoff score of 15 (all scores ≤ 10; see Table 1), and the Child Behavior Checklist for ages 1½-5 years (CBCL/1½-5; Achenbach and Ruffle, 2000), used here as a screener for broadband behavioral problems.

MRI data acquisition

All participants underwent MRI scans during natural nocturnal sleep on a GE Discovery MR750 3 T MRI scanner at the UCSD Center for Functional Magnetic Resonance Imaging, using a Nova Medical 32-channel head coil. A multiband multi-echo planar imaging (EPI) sequence allowing simultaneous acquisition of multiple slices was used to acquire two fMRI runs (400 volumes per each 6-minute run) with high spatial resolution and fast acquisition (TR = 800 ms, TE = 35 ms, flip angle = 52°, 72 slices, multiband acceleration factor = 8, 2 mm isotropic voxel size, matrix = 104 × 104, FOV = 20.8 cm). Two separate 20 s spin-echo EPI sequences with opposing phase encoding directions were also acquired using the same matrix size, FOV, and prescription to correct for susceptibility-induced distortions. High-resolution anatomical images were acquired with a fast 3D spoiled gradient recalled (FSPGR) T1-weighted sequence (0.8 mm isotropic voxel size, NEX = 1, TE/TI = min full/1060 ms, flip angle = 8°, FOV = 25.6 cm, matrix = 320 × 320, receiver bandwidth 31.25hz). Motion during anatomical scans was corrected in real-time using three navigator scans and real-time prospective motion correction (PROMO) (White et al., 2010), and images were bias corrected using the GE PURE option.

In preparation for the scan night, and to optimize MRI data acquisition, a comprehensive sleep MRI habituation protocol was implemented. An individualized scan night sleep strategy (e.g., time of arrival, approximating home-like sleeping arrangements, including access to a double MRI bed for co-sleeping families, rocking chair, modular playpen mounted on the MRI bed, lighting in the MRI suite, etc.) was developed for each child, based on the typical bedtime routines and habits assessed in advance with an in-house developed sleep habits questionnaire. To habituate the child to the scanning environment, the parents were instructed to practice inserting soft foam child-size earplugs after the child had fallen asleep, and to play an mp3 file containing the MRI sounds of the scan sequences employed in the study at progressively louder volumes for a week. On the night of the scan, noise protection was achieved with child-size earplugs and MRI compatible sound reducing headphones. Scanning commenced after approximately 30−50 min of sleep, with fMRI scans acquired first, followed by the T1 sequence.

Data analysis

MRI data pre-processing

MRI data were preprocessed with FMRIB’s Software Libraries (FSL v5.0.10) (Smith et al., 2004), Matlab 2015b (Mathworks Inc., Natick, MA, USA) using SPM12 (Wellcome Trust Centre for Neuroimaging, University College London, UK), and the CONN toolbox v17f (Whitfield-Gabrieli and Nieto-Castanon, 2012) (http://www.nitrc.org/projects/conn). Preprocessing steps included correction for susceptibility-induced distortions using the two spin-echo EPI acquisitions with opposite phase encoding directions and FSL’s TOPUP tools; motion correction using rigid-body realignment as implemented in SPM12; spatial smoothing using a 6 mm Gaussian kernel at full-width half maximum; outlier detection using the Artifact Detection Toolbox as installed with CONN v17f (ART; https://www.nitrc.org/projects/artifact_detect) with outliers defined as volumes with frame-wise displacement (FD) >0.5 mm and/or changes in signal intensity > 3 standard deviations; nuisance regression including censoring of the outliers detected by the ART toolbox and regressing the 6 motion parameters and their derivatives, and the first five PCA components derived from the CSF and white matter signal using aCompCor (Behzadi et al., 2007); and band-pass temporal filtering (0.008–0.08 Hz).

The structural images were co-registered to the mean functional image, segmented and normalized to the Montreal Neurological Institute (MNI) atlas space using non-linear registration and the default tissue probability maps included with SPM12 (for a detailed discussion on spatial normalization see Supplementary Methods). The white matter (WM) and CSF probability maps obtained from segmentation of the structural image for each subject were thresholded at 0.95 and eroded by 1 voxel. These thresholded and eroded masks were applied to functional images to extract white matter and CSF time courses, which were submitted to a principal component analysis with aCompCor (Behzadi et al., 2007) for subsequent nuisance regression. Functional images were directly normalized to MNI space with the same non-linear registration as used for the structural images.

In order to ensure that the findings were not due to individual differences in motion, mean head motion indexed by root mean square of displacement (RMSD) across two fMRI runs, was included as a covariate in all imaging analyses. RMSD calculated from rigid-body realignment of the raw data prior to TOPUP correction and the percentage of censored volumes averaged across two fMRI runs are reported in Table 1.

Independent component analysis

Preprocessed fMRI data (concatenated across two runs) from all participants were submitted to group independent component analysis (ICA) using CONN’s ICA implementation (Calhoun et al., 2001) to generate intrinsic functional networks. Each subject contributed 800 volumes (across two runs) to the group-ICA for a total of 37,600 3D volumes. Twenty independent components (ICs) were extracted, and each component’s spatial distribution and time course were visually inspected by two independent raters. ICs identified as noise (i.e., motion, cerebral spinal fluid pulsations, signal from large blood vessels) were discarded. The remaining ICs were compared to the 10 components generated by Smith et al. (2009) and the 8 components generated from the Human Connectome Project Consortium’s 500 Subjects Release (https://db.humanconnectome.org/data/projects/HCP_500), both based on adult data, as well as to published RFNs in young school-age children (Muetzel et al., 2016; Thornburgh et al., 2017). This resulted in 17 ICs being classified as RFNs (for details on the two-step IC identification and selection process, see Supplementary Methods).

Functional network metrics: within-RFN connectivity and spatial variability

Dual regression as implemented in the CONN toolbox was applied to obtain the subject-specific time courses and spatial maps corresponding to the 17 retained ICs. The group ICs of the selected networks were thresholded at z>3.0, masked with the MNI gray matter mask, and extracted as group RFN maps (see Supplementary Table S1 for coordinates and cluster description). The clusters within the thresholded RFN maps were used as regions of interest to extract the within-RFN connectivity and spatial variability measures from the corresponding subject-specific spatial maps. Within-RFN connectivity was calculated as the mean of the subject’s z-statistic for each RFN, masked by the corresponding group RFN map. This metric indicates the degree of correspondence between the individual RFN and the group RFN maps. Within-RFN spatial variability was calculated as the standard deviation of the z-statistic for each RFN, masked by the corresponding group RFN map. This measure indicates how spatially variable or homogenous the RFN is in an individual subject.

Associations with age and developmental indices

In order to examine age-related effects in the RFNs identified in ICA, multiple linear regression models were employed, with the within-RFN connectivity or spatial variability as dependent variables and age and head motion as independent variables. Correction for multiple comparisons was achieved using Benjamini-Hochberg False Discovery Rate (FDR) correction as implemented in MATLAB (corrected pFDR<0.1). Multiple linear regression analyses were also conducted to examine the relationship between indices of functional network organization and overall developmental skills indexed by the MSEL Early Learning Composite score, while controlling for age and head motion, with Benjamini-Hochberg FDR correction applied to control for multiple comparisons (corrected pFDR<0.1).

Results

Resting-state functional networks (RFNs) identified in the sample

The 17 components identified as resting-state functional networks were largely consistent with the RFNs previously described in adults and children (see Fig. 1 for spatial maps, and Supplementary Table 1 for detailed information regarding the major clusters constituting each RFN). Among the 17 RFNs, nine were unimodal primary sensory components, including sensorimotor (motor, lateral somatosensory, and medial somatosensory), auditory, and visual (medial, lateral, and occipital pole components and two higher-order visual processing) networks; six higher-order networks, including the default mode (with two subcomponents), salience, dorsal attention, right frontoparietal, and language networks; and two additional RFNs corresponding to cerebellum and subcortex (brainstem and thalamus). When compared to the previously published adult RFNs and networks obtained in school-age children (Muetzel et al., 2016; Smith et al., 2009; Thornburgh et al., 2017), the current set of RFNs notably includes distinct motor (SM1) and two somatosensory (lateral and medial; SM2, SM3) networks, which have been subsumed in one single network in past studies with school-age children and adults (e.g., Muetzel et al., 2016; Smith et al., 2009; Thornburgh et al., 2017). We also identified two DMN-like subcomponents, consistent with other reports in young school-age children (between 6 and 10 years old; Muetzel et al., 2016; Thornburgh et al., 2017), in contrast to findings in adults with the DMN typically identified as a single component with all its canonical nodes (e.g., Smith et al., 2009).

Fig. 1

Intrinsic functional networks in typically developing 1-to-3-year-old children. Results of the 20-dimensional group ICA; images are z statistics thresholded at z = 3.0 (p < .001) grouped into functional domain categories as depicted. Images are presented in the Montreal Neurological Institute (MNI) space, in neurological convention (with the left side of the brain represented on the right).

Age-related effects in network indices

Linear regression analyses of within-RFN connectivity and spatial variability with age revealed that the within-RFN spatial variability of one of the DMN components (DMN2) was negatively associated with age, after controlling for head motion and correcting for multiple comparisons (standardized coefficient = −0.59, puncorr = 0.002, corrected pFDR = 0.03), with lower within-DMN2 spatial variability observed in older children (Fig. 2). No significant associations were observed between age and other higher-order networks, for either within-RFN connectivity or spatial variability, after adjusting for head motion and multiple comparisons (see Supplementary Figure S3 for effect sizes, and Supplementary Table S2 for a full set of regression results). As predicted, no associations with age were detected for any of the primary sensory networks (see Supplementary Table S2).

Fig. 2

DMN2 iFC spatial variability across age (partial regression plot, controlling for head motion). The values on the X and Y axes reflect residuals of age in days and standard deviations of DMN2 spatial maps, after controlling for head motion (RMSD).

Within-network spatial variability and its links with developmental indices

Linear regression analyses examining links between within-RFN spatial variability and overall developmental skills indexed by the MSEL Early Learning Composite (ELC) score (corrected pFDR < 0.1) revealed that greater ELC scores were associated with lower within-RFN spatial variability of the posterior DMN component (DMN1; standardized coefficient = −0.55, puncorr = 0.008; Fig. 3A) and of the salience network (SN; standardized coefficient = −0.59, puncorr = 0.001; Fig. 3B), after controlling for age and head motion. (For a full set of multiple linear regression analyses examining links between within-RFN spatial variability in all cortical networks and MSEL ELC, see Supplementary Table S3.)

Fig. 3

Associations between iFC spatial variability in higher-order RFNs (default mode, salience, and language networks) and developmental skills. Partial regression plots between standard deviations of default mode network (A) and salience network (B) spatial maps and MSEL Early Learning Composite (ELC) scores. Panel C depicts a partial regression plot of the association between standard deviations of language network spatial maps and MSEL language index (average of MSEL Receptive and Expressive Language T scores). The values on the X and Y axes reflect residuals of MSEL scores and spatial maps standard deviations, after controlling for head motion (RMSD) and age. Greater MSEL scores correspond to more advanced developmental skills.

Finally, given the prominence of language skills in the overall cognitive and social development at this age range and because of the direct relevance of the language RFN to the emerging language skills, we examined, in a post-hoc analysis, the relationship between within-RFN spatial variability of the language network and children’s language scores measured with MSEL Verbal index (obtained by averaging Receptive and Expressive Language T scores). Greater language skills indexed by MSEL Verbal index were associated with lower within-RFN spatial variability in the language network (standardized coefficient = −0.47, p = 0.02) after controlling for age and head motion (Fig. 3C).

Discussion

In this study, we examined large-scale resting-state functional networks in a cross-sectional cohort of 1-to-3-year-old typically developing children using resting-state fMRI data acquired during natural sleep. A set of RFNs, identified through data-driven group-ICA, largely corresponded with the RFNs reported in previous studies of older children and adults. Our hypothesis-driven analyses focusing on patterns observed in higher-order v. primary sensory networks revealed two principal findings. First, analysis of age-related effects in functional network organization quantified with two network indices, within-RFN connectivity and spatial variability, revealed a negative association with age in spatial variability of one of the DMN components, with older children displaying less spatially variable, or more homogeneous DMN subcomponent consisting of the anterior and posterior midline nodes. Such age-related effect was not present in any of the primary sensory networks (e.g., visual or sensorimotor networks). Second, lower spatial variability within the posterior default mode and salience networks was associated with greater overall developmental skills, and lower spatial variability within the language network was associated with greater language skills, while controlling for age.

Negative association with age in network spatial variability: Sign of increasing specialization?

While the finding of a negative association with age in within-DMN spatial variability is consistent with our hypothesis, the predicted positive relationship between age and within-RFN connectivity in higher-order networks was not detected. Together, these results suggest that the two network indices – within-RFN connectivity and spatial variability – likely reflect distinct aspects of network organization and maturation and may unfold along a different developmental timeline in the first years of life. Within-RFN connectivity represents the extent to which the individual spatial map corresponds to, or matches the group RFN map, which can serve as a proxy for average “connectivity” within the network. On the other hand, within-RFN spatial variability, quantified as standard deviation of the z statistic, indexes the degree to which each RFN is spatially homogenous or cohesive in an individual subject, which can be conceptualized as a measure of network homogeneity or spatial specialization. Given this distinction between the two metrics, our results suggest that while the average connectivity in higher-order networks may be relatively stable across the sampled age range, between ages 1.5 and 3.5 years, these supra-modal, associative networks continue to become more spatially specialized during this pivotal developmental period. Specifically, the observed lower spatial variability (or greater homogeneity) of one of the DMN components in older toddlers indicates a shift from a more variable spatial signature to a more constrained and possibly more specialized functional organization of the DMN as children become older. This interpretation is largely consistent with the Interactive Specialization theory of functional brain development (Johnson, 2011), which proposes that brain functions – broadly construed – become increasingly specialized with development, with brain activation patterns becoming more restricted to more circumscribed brain regions (i.e., spatial specialization), as well as being elicited by a more circumscribed set of conditions (i.e., functional specialization). As cortical regions become more specialized and their response patters more finely tuned, so do the larger-scale functional networks, which represent coordinated neural activity among interconnected regions with history of common experience-driven co-activation. In support of this theory, it has been shown that, over the course of childhood, brain networks become increasingly segregated, both spatially and with regard to their functionality (Durston et al., 2006; Fair et al., 2007), reflecting growing specialization or “focalization” across age (Dehaene-Lambertz et al., 2018; Richardson et al., 2018). Thus, reduced spatial variability or increased network homogeneity observed in our sample for one of the DMN components may reflect this experience-based specialization, signifying ongoing DMN functional maturation.

The DMN is generally associated with internally oriented, self-referential mental activity, such as autobiographical memory (Svoboda et al., 2006) and mental time travel (Spreng et al., 2009), as well as social-referential cognition (Mars et al., 2012), although how (or whether) these DMN-associated cognitions manifest in young children, whose internal states we cannot directly discern, is unknown. The default mode network mainly encompasses the posterior cingulate cortex (PCC), medial prefrontal cortex (mPFC), lateral temporal cortex (LTC), and inferior parietal lobule (IPL; Buckner et al., 2008; Fox and Raichle, 2007). A primitive and fragmented DMN can already be identified in neonates and it becomes more complete and inter-connected over the course of development (Fair et al., 2008; Gao et al., 2009). For instance, Gao et al. (2009) found an increase in the number of significant connections between distinct DMN nodes from early infancy (in neonates) through the 1st year, and a subsequent decrease from age one to two years, which suggests a potential specialization process, whereby redundant functional connections are being eliminated throughout early childhood. Notably, only one of the DMN components identified in our study showed reduced within-RFN spatial variability across age. This DMN component (DMN2) includes a prominent mPFC / anterior cingulate cluster, which is not part of the other, posterior DMN component (DMN1). This finding is consistent with the more updated view of the DMN comprising at least two separate (sub-)networks with spatial division along the posterior-anterior midline (Buckner and DiNicola, 2019). Our results may thus be interpreted in support of the differential maturational timeline of the two DMN sub-networks, with the posterior DMN undergoing spatial and perhaps functional specialization ahead of the anterior DMN.

Links with developmental skills

In addition to examining age-related effects in network organization indices of RFNs in toddlerhood, another goal of our study was elucidating potential links between indices of functional network organization and emerging developmental skills. Our results revealed a consistent pattern of greater network homogeneity in higher-order RFNs associated with more advanced developmental skills, regardless of age. Namely, we found that lower within-RFN spatial variability of the posterior DMN and salience network was related to higher overall developmental skills, indexed with a standardized measure of early learning, considered a precursor of intellectual abilities (which cannot be reliably established until later in childhood; Stein and Lukasik, 2009; Bishop et al., 2011). Further, in a domain-specific comparison, lower within-RFN spatial variability of the language network was associated with higher language skills, while controlling for age. This pattern suggests that greater spatial cohesiveness (suggesting increased functional specialization) of the higher-order functional networks corresponds with more favorable developmental outcomes in toddlers. Some recent studies have shown correspondence between brain functional network connectivity and emerging developmental skills, including gross motor function and social communication skill (i.e., initiation of joint attention) in infants with high and low familial risk for autism spectrum disorder (Eggebrecht et al., 2017; Marrus et al., 2018). In a cohort of typically developing young children across a relatively wide age range (between 3 months and 6 years of age), Bruchhage et al. (2020) also found associations between functional connectivity and motor, non-verbal reasoning, and language skills. However, these previous studies used indices of network maturation that were derived based on the predefined regions of interest obtained from adult studies. Given the remarkable cognitive and functional network development between early childhood and adulthood, this strategy may have led to biased or skewed-toward-adult estimates of functional connectivity and functional networks. In contrast, our study used a data-driven approach to derive large-scale functional networks using fMRI data from a cohort of typically developing toddlers.

The relationship observed between language skills and within-RFN spatial variability of language network is particularly noteworthy because language emergence is one the most critical developmental milestones in toddlerhood. Estimated to be present in approximately 95 % of right-handed adults (Bookheimer, 2002; Szaflarski et al., 2006), left-lateralized organization of language network is thought to emerge in the first year of life, with the functional asymmetry continuing to increase with age and becoming comparable to that of adults by early school age (Berl et al., 2014; Dubois et al., 2016; Geng et al., 2012; Weiss-Croft and Baldeweg, 2015). Consistent with the notion of the ongoing lateralization (and spatial and functional refinement) of the language network in toddler years, we identified a left lateralized language network, which included canonical language regions such as planum temporale, and pars opercularis and pars triangularis of the inferior frontal gyrus, along with superior frontal regions not typically associated with language processing. Importantly, our results revealed a significant relationship between reduced spatial variability (or greater spatial cohesiveness) in this network and overall language skills, independent of age. These results suggest that the canonical language network encompassing key frontal-temporal regions is not fully integrated yet at this age, and that functional connectivity in bilateral frontal-temporal language regions may be critical for the development of language skills during toddler years.

Among the limitations of the current study is its relatively modest sample size, due to the challenges of acquiring low-motion imaging data in this age group, and the use of cross-sectional data to investigate age-related effects in brain functional organization. Critically, future longitudinal work in larger samples is needed to elucidate within-subject trajectories of functional network development and their relationship with developmental skills. Besides the effect of age on functional network organization, other factors such as socioeconomic factors and early life adversity, which are not directly assessed in the current study, may also impact the trajectory of functional network development (cf. Olson et al., 2021).

In conclusion, our results provide the important initial evidence of age-related effects in brain functional organization in a cross-sectional cohort of typically developing toddlers. The current study also identified robust relationships between greater network spatial cohesiveness, thought to reflect increasing functional specialization, and more advanced developmental skills.

Data availability

All behavioral and imaging data from this study will be shared through the National Institute of Mental Health Data Archive (NDA), Study Collection #2338. The final ICA group components (RFN maps) generated in this cohort of typically developing young children are available for download on Open Science Framework (https://osf.io/nkmw9/; DOI 10.17605/OSF.IO/NKMW9).

Declaration of Competing Interest

All the authors declare no conflicts of interest.