The emergence of a functionally flexible brain during early infancy.

Adult brains are functionally flexible, a unique characteristic that is thought to contribute to cognitive flexibility. While tools to assess cognitive flexibility during early infancy are lacking, we aimed to assess the spatiotemporal developmental features of "neural flexibility" during the first 2 y of life. Fifty-two typically developing children 0 to 2 y old were longitudinally imaged up to seven times during natural sleep using resting-state functional MRI. Using a sliding window approach, MR-derived neural flexibility, a quantitative measure of the frequency at which brain regions change their allegiance from one functional module to another during a given time period, was used to evaluate the temporal emergence of neural flexibility during early infancy. Results showed that neural flexibility of whole brain, motor, and high-order brain functional networks/regions increased significantly with age, while visual regions exhibited a temporally stable pattern, suggesting spatially and temporally nonuniform developmental features of neural flexibility. Additionally, the neural flexibility of the primary visual network at 3 mo of age was significantly and negatively associated with cognitive ability evaluated at 5/6 y of age. The "flexible club," comprising brain regions with neural flexibility significantly higher than whole-brain neural flexibility, were consistent with brain regions known to govern cognitive flexibility in adults and exhibited unique characteristics when compared to the functional hub and diverse club regions. Thus, MR-derived neural flexibility has the potential to reveal the underlying neural substrates for developing a cognitively flexible brain during early infancy.

Resting-state functional MRI (rsfMRI), an imaging method characterizing temporal synchrony of blood oxygen level-dependent (BOLD) contrast signals among different brain regions, has been widely employed to delineate brain functional network topologies. Because rsfMRI acquires time-series images when subjects are in a resting condition without any goal-directed task demands, it is an optimal tool to reveal brain functional networks in pediatric subjects who are unable to comply with the requirements of performing tasks for goal-directed fMRI. Extensive results using rsfMRI have been reported in early brain development. The developmental pattern of the widely reported functional networks in adults follows a temporal sequence from maturation of the primary sensory networks to higher-order functional networks (1). The emergence of the default mode network prenatally (2) and postnatally (3), the presence of a small world topology of the brain networks in infants (4), and three distinct developmental stages (i.e., 0 to 1 mo, 2 to 7 mo, and 8 to 24 mo) (5) have been observed. In addition, the number of functional modules increases throughout the first year of life, with strengthened intra- and intermodular connections, thought to support functional segregation and integration (6). The emergence of functional hub regions in the brain (i.e., regions critical for efficient communication) has also been reported from birth to 2 y of age (4). While the aforementioned studies have offered invaluable insights into early brain functional developmental processes, they assume that functional network organization is static across the course of a resting-state scan. However, there is increasing evidence in adults that functional network organization is time-varying, and that dynamic network processes underlie meaningful aspects of cognition and behavior (7). A key finding from these studies is that functional brain network organization dynamically changes with learning. For example, using a sliding window approach to analyze rsfMRI time-series images, Bassett et al. (8) quantified dynamic modular reconfiguration when adult participants were learning a new motor skill. Specifically, they derived a normalized parameter they termed “functional flexibility,” which quantified the frequency at which nodes (brain regions) changed their allegiance from one functional module to another. They reported that functional flexibility first increased and then decreased while learning a new motor skill. In addition, functional flexibility in one learning session predicted the amount of learning in a future session.

In the context of early brain development, measures of functional flexibility during the first years of life could offer the opportunity to assess reorganization of brain functional network topologies that may result from learning via interaction with the external environment. In particular, a plethora of studies have demonstrated that external stimuli and learning shape the maturation processes of early brain functional development (9–13). Moreover, the definition of functional flexibility, which assesses the frequency at which a brain region changes its functional role (i.e., changes its allegiance from one functional module to another), resembles the classic definition of cognitive flexibility, or the readiness with which one can selectively switch between mental processes to appropriately respond to environmental stimuli (14). Cognitive flexibility, one of the three main components of the executive functions (15, 16), is a critical feature of human cognition (17) and has been reported to predict reading ability (18) and future academic success (19). In addition, higher cognitive flexibility is associated with higher resilience to stress and better creativity in adulthood (20–22). In contrast, poor cognitive flexibility has been used as a biomarker for brain diseases, including autism (17), Parkinson’s disease (23), Alzheimer’s disease (24), schizophrenia, and mania (25). Recent studies suggest that brain variability is an important feature underlying cognitive flexibility (22, 26, 27); thus, understanding how variability, or flexibility in network organization, matures across early development could elucidate the development of this critical cognitive process. To better reflect the neural basis of functional flexibility and distinguish it from cognitive flexibility, “neural flexibility” is used hereafter to represent rsfMRI-derived functional flexibility. The ability to noninvasively measure neural flexibility is likely to open a window of opportunity to assess developmental features and the emergence of a functionally flexible brain during early brain development, an area of research that could yield highly profound insights into the development of higher-order brain functions and that could ultimately be used to predict cognitive outcome.

We leveraged a longitudinal imaging study in typically developing children, the Multivisit Advanced Pediatric Brain Imaging Study for Characterizing Structural and Functional Development (MAP Brain Imaging Study) (1, 5, 6, 28), to uncover the developmental features of neural flexibility at both the regional and network levels to characterize the emergence of a functionally flexible brain during early infancy. Since increased neural flexibility has been reported to be associated with learning new skills in adults (8), and learning/interaction with the external environment plays a vital role in early brain development, we hypothesized that spatial and temporal heterogeneity of neural flexibility may represent varying maturation patterns of different brain functional domains. We further hypothesized that the neural flexibility trajectories of higher-order brain functional networks would exhibit a faster pace of increasing neural flexibility than that of basic brain functional networks. Finally, we hypothesized that neural flexibility in early infancy could predict infants’ behavioral ability at a later age. This study further reported how quantitative measures of neural flexibility could shed light on the development of cognitive flexibility through defining the “flexible club” as brain regions exhibiting significantly higher neural flexibility than that of the whole brain. These regions therefore represent brain regions that may be critical for state-related transitions, including transitions between cognitive tasks with different demands (i.e., cognitive flexibility).

Results

Typically developing children (n = 52) were densely and longitudinally MR-imaged (up to seven times) during natural sleep without sedation over the first 2 y of life (Fig. 1 ). The number of subjects and sex information for each age group is provided in Table 1. No significant difference was found for the framewise displacement (FD) of rsfMRI data among age groups [ANOVA F test, F(6,196) = 1.8, ]. In addition to undergoing longitudinal MR imaging sessions, the General Conceptual Ability (GCA) score of the Differential Ability Scales II (DAS) test was obtained from 31 subjects at 5/6 y of age. The numbers of subjects whose DAS scores were available for each age group are also provided in Table 1. Using a sliding window approach, a multilayer network structure was constructed over all time windows and the GenLouvain (29) algorithm was used to evaluate dynamic network reconfigurations (Fig. 1). Subsequently, nodal temporal allegiance to functional modules was calculated (Fig. 1), which in turn provided quantitative measures of regional neural flexibility.

Fig. 1.

Illustration of data distribution and multilayer community detection. (A) The distribution of the age information for all subjects. (B) Examples of T1-weighted and T2-weighted images obtained from a child who underwent longitudinal MRI at birth and 3, 6, 9, 12, 18, and 24 mo of age. (C) Illustration of the multiplayer networks. Each node was linked to itself one time window before and after. (D) Examples of dynamic module transition for nodes with frequent module transitions (Top) and few transitions (Bottom).

Table 1.

Number of subjects and motion parameters of each age point

	Age, mo
	<1	3	6	9	12	18	24
No. of subjects	29	26	33	32	32	31	20
Male	15	14	17	17	16	17	8
Female	14	12	16	15	16	14	12
FD (mean ± SE), mm	0.22 ± 0.04	0.16 ± 0.03	0.15 ± 0.02	0.12 ± 0.01	0.15 ± 0.03	0.21 ± 0.05	0.09±0.01
DAS evaluated at 5/6 y	16	16	17	19	19	22	17

Spatiotemporal Characteristics of Brain Neural Flexibility.

As indicated above, increased neural flexibility has been associated with learning new skills and predicting learning ability in adults (8). In the context of early brain development when learning and interaction with the environment play a vital role, we tested whether brain neural flexibility increased with age. While a substantial spatial variability of neural flexibility is present (Fig. 2), significantly elevated global neural flexibility was found from birth to 2 y of age [linear mixed-effects (LME) model, F test, F(1,201) = 25.13, ; Fig. 2]. Due to the reported difference in maturation trajectories of different functional systems (30, 31), we tested whether there was heterogeneity among the development patterns of neural flexibility in different regions and systems. The brain regional slopes of neural flexibility trajectories are shown in Fig. 2. Many of the regions in the frontal lobe exhibited a higher slope than other brain regions. To further probe the neural flexibility patterns in the core regions of several commonly reported brain functional networks, the developmental patterns from two primary brain functional regions (calcarine fissure for visual and precentral gyrus for motor) and seven higher-order brain regions (precuneus cortex, inferior temporal gyrus [ITG], anterior cingulate cortex [ACC], inferior frontal gyrus [IFG], temporoparietal junction [TPJ], posterior parietal cortex [PPC], and anterior insula [aI]) over the first 2 y of life were chosen and are shown in Fig. 2. The anatomical locations of these regions are shown in Fig. 2 (white filled circles). These higher-order brain regions were chosen because they have been implicated in cognitive flexibility (22) (ACC, IFG, TPJ, PPC, and aI), the default mode network (precuneus) (32) and higher-order visual function (ITG) (33).

Fig. 2.

Spatiotemporal development of brain neural flexibility. (A) Spatiotemporal distribution of regional neural flexibility at 0, 3, 6, 9, 12, 18, and 24 mo of age. (B) The developmental pattern of whole-brain neural flexibility. The red line and shaded area represent the fitted trajectory of whole-brain neural flexibility and 95% confidence interval (LME model, F test, ). The blue boxplots represent the distribution of whole-brain neural flexibility at each age group. (C) The developmental slope of the neural flexibility at each brain region from birth to 2 y of age is shown on brain surfaces (colored bar represents the slope values). (D). Regional-specific developmental patterns of brain regions, including calcarine fissure, precentral gyrus, precuneus, ITG, ACC, IFG, TPJ, PPC, and aI. These chosen brain regions are marked with white filled circles in C. Blue boxplots represent the distribution of neural flexibility of each region at 0, 3, 6, 9, 12, 18, and 24 mo of age. The P values of the LME fitted trajectories after FDR correction were as follows: , , , , , , , , and .

Distinct temporal patterns were observed between the calcarine fissure (primary visual cortex) and precentral gyrus (primary motor cortex). The calcarine fissure exhibited a slight decrease [LME model, F test, F(1,201) = 3.63, , false discovery rate (FDR)-corrected], whereas the precentral gyrus exhibited an increase [LME model, F test, F(1,201) = 19.82, , FDR-corrected] in neural flexibility with age. In contrast, the temporal patterns of higher-order brain functional areas were more variable. The precuneus cortex had a marked decrease from birth to 3 mo [two-sample t test, t (53) = −2.8, , uncorrected], followed by a steady increase of neural flexibility with age [LME model, F test, F(1,201) = 10.49, , FDR-corrected]. The ITG showed persistently high neural flexibility across all ages without an age effect. For cognitive flexibility-associated regions, neural flexibility increased substantially in the ACC [LME model, F test, F(1,201) = 10.93, FDR-corrected], TPJ [LME model, F test, F(1,201) = 9.15, FDR-corrected], and aI [LME model, F test, F(1,201) = 8.29, , FDR-corrected], while a gradual increase, although not significant, was found in the IFG [LME model, F test, F(1,201) = 4.08, , FDR-corrected] and PPC [LME model, F test, F(1,201) = 2.72, , FDR-corrected]. Additional results on neural flexibility changes between two adjacent age groups are provided in .

We further evaluated network-level developmental trajectories, including those of representative unimodal functional networks (the visual and motor networks) and those implicated in cognitive flexibility in adults (22) (the salience network, dorsal attention network, ventral attention network, working memory network, and inhibition network) using the Shen268 functional parcellation (34) (Fig. 3). More information on the use the Shen268 atlas and the selection of these networks is provided in . Additional results using other higher-order brain functional networks are shown in . The anatomical locations of these networks are summarized in Fig. 3 and .

Fig. 3.

Development of brain neural flexibility at network level and the prediction of brain neural flexibility to later behavioral performance. (A) Developmental patterns of neural flexibility in visual and motor-related networks (blue shaded area) and networks associated with cognitive flexibility (orange shaded area). The anatomical features of these brain networks are shown (Upper). For each network, the red line and shaded area represent the fitted trajectory of neural flexibility and 95% confidence interval, respectively. The P value of the LME-fitted trajectories after FDR correction were as follows: , , , , , , , , and . The blue boxplots represent the distribution of neural flexibility of each brain network at each age group. DAN: dorsal attention network; VAN: ventral attention network; WM: working memory network. (B) The developmental slopes of the neural flexibility of whole brain and different brain networks, sorted by the slope values. The $, #, *, ^, and & symbols indicate significantly different in slopes (LME model, F test, P < 0.05, FDR-corrected) with the visual I, visual association, whole brain, motor, and inhibition networks, respectively. (C) The prediction of brain neural flexibility measured during early infancy to later behavioral performance. Neural flexibility of visual I network at month 3 was significantly correlated with GCA score evaluated at 5/6 y of age [Pearson’s correlation t test, t (14) = −5.15, , FDR-corrected].

It has been widely documented that primary visual function matures early in infancy (1) to enable the visual processing capability observed in neonates (35) and that these regions are mostly associated with unimodal functions (36). Consistent with these findings, we observed that the neural flexibility of visual-related networks exhibited no age effects, whereas significant age effects were observed for the developmental trajectories of all remaining networks [Fig. 3; LME model, F test, Fmotor(1,201) = 12.18, Fsalience(1,201) = 34.41, FDAN(1,201) = 16.26, FVAN(1,201) = 9.3, FWM(1,201) = 32.67, Finhibition(1,201) = 7.89, , FDR-corrected]. In addition, the primary visual network showed a significantly lower slope than that of the dorsal attention, working memory and salience networks [LME model, F test, FV1_DAN (1,201) = 11.74, FV1_WM (1,201) = 18.05, FV1_Salience (1,201) = 23.26, , FDR-corrected] (Fig. 3). Together, these results suggest that brain regions that mature early may also reach mature levels of neural flexibility earlier, thus exhibiting a reduced change in neural flexibility over the first 2 y of life as compared to functional networks that undergo greater maturation.

Although the motor network is typically considered a primary functional network similar to that of the primary visual network, the developmental trajectory of the motor network was different from that of the visual networks. The primary motor network started with an average neural flexibility immediately after birth, followed by a significant increase in neural flexibility during the first 2 y of life (Fig. 3) [LME model, F test, F(1,201) = 12.18, , FDR-corrected]. As shown by Bassett et al. (8), increased neural flexibility is associated with greater improvements in learning a new motor skill. Thus, this continuing increase of motor network neural flexibility may be indicative of continuing development of motor skills during early infancy. Motor developmental milestones during the first years of life have been widely documented, including the development of gross motor ability of lying movement (months 2 to 4), sitting (month 6), standing (month 9), walking (month 12), and running (month 24). Fine motor ability also develops extensively in the first 2 y of life, including holding and shaking toys (month 4), as well as drinking from a cup and eating with a spoon (month 18). Thus, these results indicate that increases in motor neural flexibility may underlie the continuing development of motor skills.

The brain networks associated with cognitive flexibility and other representative higher-order brain functional networks (Fig. 3 and ) all exhibited a significant increase in neural flexibility with age [LME model, F test, Fsalience(1,201) = 34.41, FDAN(1,201) = 16.26, FVAN(1,201) = 9.3, FWM(1,201) = 32.67, Finhibition(1,201) = 7.89, Fmedialfrontal(1,201) = 22.46, Ffrontoparietal(1,201) = 19.03, Fdefault(1,201) = 20.64, , FDR-corrected]. Furthermore, their slopes were faster than that of the whole brain (Fig. 3 and ), particularly the salience network [LME model, F test, F(1,201) = 14.95, , FDR-corrected]. Notably, cognitive flexibility-related networks displayed the highest slope, suggesting the potential emergence of rudimentary cognitive flexibility during early infancy and underscoring the potential significance of assessing neural flexibility during early infancy.

Finally, a direct comparison of the trajectory of each network derived from neural flexibility to trajectories derived from functional connectivity strength, a widely used approach to characterize early brain functional development (1, 28, 37), is provided in . This comparison suggests that the trajectories of brain neural flexibility may be more representative of brain functional maturation during early brain development.

Early Brain Neural Flexibility Is Associated with Behavioral Assessments at Later Ages.

A large intersubject variability of neural flexibility was noted within each age group (Figs. 2 and 3), which suggests differences of maturation rates among subjects. Since neural flexibility has been reported to predict learning outcome in adults, we tested whether neural flexibility in early infancy could predict infants’ behavioral ability at a later age. In addition to undergoing longitudinal MRI sessions, the GCA score of the DAS test was obtained at 5/6 y of age in 31 subjects. Correlation analyses were performed between network-level neural flexibility and GCA. A negative correlation was observed between visual I neural flexibility at 3 (Fig. 3) and 18 mo () of age with later GCA [Pearson’s correlation t test, t03mon (14) = −5.15, t18mon (20) = −2.08, , ]. After multiple comparison correction, the visual I network at month 3 remained significantly correlated with GCA score (, FDR-corrected). No relationships between neural flexibility of other networks and GCA were significant. Furthermore, consistent with the network-level results, regional analyses revealed that neural flexibility of multiple regions in the primary visual areas was negatively associated with GCA at several age points, whereas some regions associated with cognitive flexibility were positively associated with GCA scores () (Pearson’s correlation t test, , uncorrected).

Flexible Club.

Regions with high neural flexibility actively change their affiliation among functional modules. Therefore, they may provide the ability to quickly change functional organization in response to external stimuli. There is evidence that multimodal brain functional areas have a high level of neural flexibility (38), potentially facilitating task transitions over time. Thus, this set of regions could be critically important for cognitive flexibility (22, 39, 40).

To identify brain regions with high neural flexibility, we defined the “flexible club” as the group of brain regions with neural flexibility that is significantly higher than that of the whole brain (Fig. 4) (details of the definition of flexible club are provided in ). We hypothesized that the brain flexible club would include brain regions essential for cognitive flexibility. To delineate potential differences in the flexible club in early infancy as compared to those in adolescents and adults, flexible clubs were identified using the Adolescent Brain Cognitive Development (ABCD) (41) and the Human Connectome Project (HCP) datasets (42) (see Fig. 4 for comparisons). The flexible clubs of the ABCD and HCP cohorts comprised brain regions critical for cognitive flexibility, including anterior insula, TPJ, dorsolateral/ventrolateral prefrontal cortex, IFG, frontal eye field, and subcortical regions (Fig. 4), supporting our hypothesis that brain flexible clubs are critical for cognitive flexibility (details of the ABCD and HCP data are provided in ). In our study cohort, while the memberships of the flexible club varied with age during early infancy, several patterns were observed. First, starting from birth, anterior ITG, IFG, orbital frontal cortex, insula, and TPJ were consistently present in the flexible club in most of the age groups. Second, most of the regions observed in the adolescent and adult flexible clubs were included in the flexible club at 2 y, including ventrolateral prefrontal cortex, dorsolateral prefrontal cortex, premotor area, aI, and TPJ. As mentioned above, these regions are core regions underlying cognitive flexibility (22). Although these results cannot conclusively determine that cognitive flexibility ability is present at 2 y of age, nor do they equate the flexible club with cognitive flexibility, they indicate that regions critical for cognitive flexibility in adults are among those that have statistically higher neural flexibility in the brain during the first 2 y of life.

Fig. 4.

Brain flexible club. (A) An example of ranked regional flexibility from month 0 (red/blue indicate regional neural flexibility significantly higher/lower that that of whole-brain neural flexibility; orange indicates regions with no significant difference, paired t test, ). Error bars indicate the SEs. (B) Spatial distribution of the brain flexible clubs. The flexible clubs of adolescents (ABCD) and adults (HCP) were included for comparison. Color bars represent neural flexibility.

Flexible Club, Functional Hub, and Diverse Club.

Using graph theory and network science, brain regions that play distinct roles in brain networks have been reported. Specifically, brain functional hub regions represent brain regions with high betweenness centrality (BC), critical for efficient communication across the whole brain (4). The notion of the brain diverse club has been recently proposed, reflecting brain regions playing a critical role in functional integration (43, 44). Given the reported different functional roles, it is highly plausible that the flexible club possesses unique characteristics when compared to brain regions in the functional hub and diverse club. Here, we examined several important graph-related characteristics among brain regions that are functional hubs, as well as those in the diverse club and flexible club, including anatomical locations, topological roles, connection strengths, and temporal stability. To make comparisons between groups, the number of functional hub and diverse club regions were set to be identical to that of the flexible club at each age group.

The topological locations of the flexible club regions were mostly located in peripheral areas of communities, which is different from functional hub and diverse club regions, which were located in the topological centers of the graph (Fig. 5). Starting from birth, the membership of the flexible club largely did not overlap (59.64%) with functional hub and diverse club regions across all age groups (Fig. 5 and ). The overlap between the flexible club and functional hub regions was 11.31% and between the flexible club and diverse club regions was 18.93%. Only 10.1% of regions were associated with all three classes. These findings are consistent across different age groups, as well as different choices of sparsity (i.e., 5%, 10%, 15%, and 20%) of static connectivity matrices (). These results demonstrate that members of the flexible club are spatially distinct from functional hub and diverse club regions.

Fig. 5.

Comparisons among brain flexible club, diverse club, and functional hub. (A) Visualization of brain communities, with the corresponding neural flexibility, participation coefficient, and BC (red/blue color corresponding to high/low parametric values). (B) The spatial overlapping ratios among brain flexible club, functional hub, and diverse club at birth. The functional hub and diverse club were generated from 10% connectivity matrices. These findings were consistent across different age groups, as well as different choices of sparsity (i.e., 5%, 10%, 15%, and 20%) of static connectivity matrices (). Statistical comparisons of brain resilience to targeted attacks (C), mean dynamic connectivity strength (D), and mean temporal instability (E) among brain flexible club, diverse club, and functional hub. The results at birth are shown and additional results across different age groups, as well as different choices of sparsity (i.e., 5%, 10%, 15%, and 20%) of static connectivity matrices are shown in . Group comparisons were performed for flexible club, diverse club and functional hub. Statistical significance levels: , , .

We further employed a simulated lesioning approach to determine resilience of functional hub, diverse club, and flexible club regions, respectively. Specifically, regions in each category were removed, and the impacts to whole-brain global efficiency were estimated after each removal. Consistent with the previously reported results (4), removing functional hub regions led to a significant reduction in global efficiency (Fig. 5) [paired t test, t (28) = −28.8, , FDR-corrected]. Similar findings were observed when removing diverse club regions [paired t test, t (28) = −11.21, , FDR-corrected]. In contrast, global efficiency was not strongly affected after removing flexible club regions [paired t test, t (28) = 1.26, , FDR-corrected], suggesting that the flexible club may be more resilient to brain injury. Similar results are also observed across different age groups, as well as different choices of sparsity (i.e., 5%, 10%, 15%, and 20%) of static connectivity matrices ().

In addition, brain regions in the flexible club exhibited significantly lower dynamic connectivity strength (Fig. 5) [paired t test; flexible club vs. functional hub: t (28) = −14.12, , FDR-corrected; flexible club vs. diverse club: t (28) = −9.41, , FDR-corrected] and higher temporal instability (Fig. 5) [paired t test; flexible club vs. functional hub: t (28) = 6.57, , FDR-corrected; flexible club vs. diverse club: t (28) = 5.19, , FDR-corrected] when compared to functional hub and diverse club regions. These findings are also consistent across different age groups and different choices of sparsity (i.e., 5%, 10%, 15%, and 20%) of static connectivity matrices (). Although many factors may explain the observed differences in dynamic connectivity strength and temporal instability among functional hub, diverse club, and flexible club regions, it is plausible that the observed weaker and temporally unstable connections between brain regions in the flexible club facilitate rapid transitions of their allegiances to different functional modules.

Discussion

Our study aimed to reveal spatiotemporal processes of the development of a functionally flexible brain in typically developing children during the first 2 y of life. Our results are in support of a modular and dynamic brain network organization immediately after birth (3–6, 37). We found that whole-brain neural flexibility significantly increased with age [LME model, F test, F(1,201) = 25.13, ], demonstrating the emergence of more flexible brain function during the first 2 y of life. Regional variability of brain functional trajectories was also apparent. Specifically, higher-order brain functional areas largely exhibited a continuing increase of flexibility with age at the regional and network levels (LME model, F test, , FDR-corrected). In contrast, distinctly different trajectories between primary visual and motor functional areas (both at the regional and at the network levels) were observed. Similar to the higher-order brain functional areas, the primary motor network exhibited an increase in flexibility with age [LME model, F test, F(1,201) = 12.18, , FDR-corrected]. In contrast, the primary visual areas showed a relatively low neural flexibility that remained low throughout the first 2 y of life. It has been widely recognized that the primary visual networks mature early in infancy (1, 45). Our finding that the flexibility of the primary visual network did not change with age is consistent with already mature network functioning. Importantly, we found a significantly negative association between neural flexibility of the primary visual network at 3 mo and cognitive ability at 5/6 y of age [Pearson’s correlation t test, t (14) = −5.13, , FDR-corrected]. Thus, low neural flexibility (more matured) of the primary visual network in early infancy was associated with better cognitive outcome at 5/6 y of age. We additionally found that members of the flexible club during early infancy included brain regions associated with cognitive flexibility in adults. Finally, the brain flexible club exhibited distinctly different characteristics from brain functional hub regions and the diverse club, including different anatomical locations, topological properties, connection strengths, and temporal instability, suggesting that the flexible club possesses a unique role in brain function.

Neural Flexibility and Behavioral Outcomes.

Although our study could not reveal and did not focus on the origins of neural flexibility (), the neural flexibility of the primary visual network at 3 mo of age was negatively associated with GCA score evaluated at 5/6 y of age. As discussed previously, a low and stable neural flexibility in the visual network represents greater maturity. Therefore, our results suggest that a more mature primary visual network (i.e., lower neural flexibility) at 3 mo of age may be beneficial for long-term cognitive development. The importance of visual functions in cognitive development has been widely reported; infant visual performance, such as attention and fixation, has been shown to predict later neurocognitive development (46–48). Specifically, visual fixation in newborns is significantly correlated with visual-motor performance at 2 and 5 y of age and associated with visual reasoning at 5 y of age (47). Infant fixation duration is also associated with childhood effortful control, surgency, and hyperactivity–inattention (48). Our results thus echo results reported in previous behavioral studies concluding that visual function in early infancy is critical for future development of higher-order cognition.

Finally, in exploratory analyses we found that regions in the flexible club such as the TPJ, ITG, and superior temporal gyrus were correlated with GCA score assessed at 5/6 y of age. We cannot make too much out of these findings, since after correcting for multiple comparisons the observed associations were no longer significant. Therefore, studies with a larger sample size will be needed to determine whether these findings are significant. In summary, while more studies are needed to prospectively evaluate how neural flexibility measures during early infancy may predict later cognitive outcomes, our results reveal a potential association between neural flexibility during early infancy and later cognition.

Neural Flexibility and Cognitive Flexibility.

Cognitive flexibility, one of the three components of executive function (15, 16), has been reported to develop starting from 3 y of age with a sharp increase between 7 and 9 y of age (22, 49). However, we argue that the increased neural flexibility that we observed from 0 to 2 y of age in regions associated with cognitive flexibility and other higher-order cognitive functions may indicate early developmental processes that support the later emergence of cognitive flexibility. In particular, it has been suggested that cognitive flexibility is accomplished through multiple networks working in synergy to accomplish mental switching. Specifically, prior to switching mental functions, the salience and attention network would be engaged for awareness of and attention toward external stimuli, the working memory network would be engaged to remember the current task, and finally the inhibition network would be engaged to inhibit the current task. Inhibition develops as early as 12 mo and is largely mature by 10 to 12 y of age (22). Working memory also emerges during the first year of life and continues to improve into adolescence (50). The attention system is active from birth as a means for directing engagement with the environment and subserving learning (50). These findings suggest the possibility that cognitive flexibility may start to develop during the first years of life as well. However, the lack of tools capable of assessing these aforementioned functions during early infancy has substantially hampered our ability to gain insights into the development of cognitive flexibility. To this end, the measure of neural flexibility could be a powerful parameter enabling the assessment of brain cognitive flexibility during early infancy for the following reasons. First, recent adult task- and resting-state fMRI studies have demonstrated that brain neural flexibility is cognitively beneficial for working-memory performance and could help to predict task performance including working-memory performance (51, 52), which is one of the three components of executive functioning (15, 16, 22). Second, neural flexibility is calculated as nodal modular transition frequency among different functional modules across time, which resembles the classic definition of cognitive flexibility. Third, unlike behavioral assessments of cognitive flexibility, which are impractical during early infancy, rsfMRI can be easily applied to infants. Fourth, members of the flexible clubs in our cohort consist of brain regions (TPJ, dorsolateral/ventrolateral prefrontal cortex, aI, premotor area, and subcortical areas) that are widely implicated in cognitive flexibility in adults. Finally, the developmental trajectories of neural flexibility of higher-order functional networks (salience, attention, working memory, and inhibition networks) implicated in cognitive flexibility all show increased neural flexibility with age, indicating continuing development of these functional networks.

Importantly, the optimal neural flexibility in a healthy brain highly depends on the functional domain. It has been found that sensorimotor and visual regions have lower flexibility, while regions associated with multimodal higher-order functions have higher flexibility (38). Indeed, in our study visual regions were observed to have low neural flexibility beginning from 3 mo (), and subjects with lower flexibility in the visual network had better cognitive outcomes at 5/6 y as assessed with GCA. While motor regions showed significantly increasing flexibility during infancy, they have been reported to be associated with low flexibility in adults (38). This observed developmental feature of neural flexibility in motor regions may underlie the learning of fine motor skills during early infancy (8). Finally, the regions with the highest flexibility (i.e., flexible club) were mostly located in multimodal association regions involved in higher-order cognition generally and cognitive flexibility specifically.

Brain functional hub regions play a critical role in whole-brain information transmission. Therefore, strong and stable connections are required for these regions to ensure the stability of the whole-brain system. Injury to these regions greatly influences whole-brain global efficiency and leads to major deficits in cognitive performance (53). Moreover, since the brain diverse club is responsible for functional integration, injury to regions of the diverse club would also greatly influence whole-brain efficiency, leading to an inability to integrate information and cause cognitive impairment (54).

In marked contrast to functional hub and diverse club regions, our results demonstrated that regions high in neural flexibility—termed the flexible club—have unique properties as compared to brain functional hub and diverse club regions. We found that, across all ages, few regions in the flexible club overlapped with functional hub or diverse club regions. Compared with functional hub and the diverse club regions, flexible club regions possessed significantly weaker connection strength (paired t test, , FDR-corrected) and higher variability of connection strength (paired t test, , FDR-corrected). These results are not surprising, since these three types of regions are responsible for different functions. More importantly, recent studies have indicated that weak connections enable the system to function within many difficult-to-reach states, reflecting a capacity to adapt to novel situations by engaging mechanisms for flexible behavior (55, 56), which echoes our findings. Finally, using a simulation, we found that injury to flexible club regions did not significantly influence whole-brain global efficiency, as these regions tend to be located at the periphery of different modules and thus are not as important for global information transmission. Therefore, the flexible club may be resilient to brain injury.

Limitations.

Our study has several limitations. First, changes of brain tissue properties during early infancy may affect the accuracy of registration. Nevertheless, our team has developed image analysis tools (57) that minimize the potential inaccuracy of registration. Second, a constant echo time (TE) was used for the echo-planar imaging sequence across our entire study. Since T2* of brain tissues is longer during early infancy and progressively shortens with age, using a constant TE over all ages could impact the sensitivity of BOLD effects. While using age-specific TEs could potentially mitigate this limitation, the changes of T2* are spatially nonuniform, making it difficult to select a single TE for the whole brain. In addition, adjusting TE based on age could introduce additional biases. Third, the number of scanning volumes (i.e., 150 volumes, 5 min) is less than that suggested by the Human Connectome Protocol (42, 58). While a longer acquisition time and greater number of volumes of rsfMRI data are clearly desirable, the choice of duration reflects a compromise in order to maximize the success rate of obtaining usable resting fMRI in imaging typically developing children without sedation. Fourth, sleep information was not available, making it difficult to account for the effects of different sleep stages on functional connectivity. Nevertheless, with the approaches and results reported in this study as well as in the literature (59, 60), we believe that the effects should be minimal. More detailed discussion of these limitations are provided in .

Materials and Methods

A total of 52 typically developing infants were included in this study. Subjects were imaged at <1, 3, 6, 9, 12, 18, and 24 mo of age. Since all 52 subjects were typically developing children, no sedation was employed for imaging and all subjects were imaged during natural sleep. The experimental protocols were approved by the Internal Review Board, University of North Carolina at Chapel Hill. Written informed consent was obtained from the parents of all participants. All subjects were scanned with 3-T MR scanner (Siemens Medical Systems), to acquire T1w , T2w , and rsfMRI images (repetition time = 2 s, TE = 32 ms, 150 volumes, ). All resting functional MR images were preprocessed using the FSL (FMRIB Software Library) (61–63) and MATLAB, following the published wavelet-based preprocessing pipeline (64). After excluding the data of failed scans with excessive motion (i.e., spike percentage >5%), a total of 203 scans were included in data analyses.

To mitigate the effects of age-dependent gray matter contrast on the accuracy of registration, a pediatric-specific tissue segmentation approach utilizing both T1w and T2w images was employed to obtain tissue segmentation results (65). A longitudinal registration pipeline was used to generate the deformation field from each subject to the standard template as well as the reverse deformation field (5, 6, 66, 67). Shen268 (34) atlas was then deformed back to subject space and mean time series of each region was estimated by averaging voxel time series in each region. Additional results using AAL (68) and CC200 (69) are provided in .

A sliding windows approach was employed with a window width of 30 volumes and a step size of 1 volume, yielding a total of 111 sliding windows. The effects of window length were evaluated and results are provided in ). Interregional connection of each sliding window was captured using Pearson’s correlation. In order to remove weak and random connections, a P value for each correlation coefficient was estimated using the MATLAB function corrcoef; only connections significantly different from zero were retained (Pearson’s correlation t test, , FDR-corrected) and converted to an absolute value. Additional results when signs of connections were considered are provided in ). Overall, although some discrepancies on the spatial distributions of neural flexibility were observed using absolute, positive, and signed connectivity values, the developmental patterns of neural flexibility remain consistent.

A multilayer network structure was constructed for each scan by linking networks between sliding windows (29) (Fig. 1). Dynamic community detection was performed using the generalized Louvain method (GenLouvain) (29, 70), which includes the ability to study the dynamic community structure in a time-dependent, multilayer, and/or multiplex network. The multilayer modularity quality function (Q) is defined as

where the matrix of layer has components . Note the layer here refers to an adjacency matrix of a given sliding window. The resolution parameter of layer is ; and are the community assignments of node in layer and node in layer , respectively; is the interlayer coupling strength parameter connecting node in layers and ; is the total edge weight in the network, calculated as ; is the intralayer strength of node in layer ; the interlayer strength of node in layer is ; the strength of node in layer is ; and is the total edge weight in layer , defined as . We kept resolution parameters the same across layers and gave all interlayer connections that were present (between the same node in layers that neighbor in time) the same weight . We optimized multiplayer modularity by running GenLouvain for 100 times for each scan, since this procedure is not deterministic and each run may yield slightly different partitions of the community. In this study, resolution and coupling parameters were set to be unity, which was used in previous studies (8, 38, 51). The effects of using different resolution and coupling parameters were evaluated in ) with a wide range from 0.2 to 2, respectively.

Given the dynamic community detection results, we defined the neural flexibility of a node as the number of times that a node changed its modular assignment across the sliding windows, normalized by the total possible changes.

The developmental patterns of neural flexibility were further evaluated by an LME model and a generalized addictive mixed model (GAMM), separately (). The generalized cross-validation (GCV) errors were used to compare the fitting performance between the two models. The differences of GCV were less than 5% in all network and 10% for all region-of-interest comparisons, suggesting that the LME and GAMM models performed similarly. Since the LME model could provide quantitative measure such as the slope to represent developmental pace while the GAMM model could not, we have chosen to report the LME model in our study. We further determined if the slopes across different functional systems were different by fitting their neural flexibility differences with age using an LME model.

In addition, the spatiotemporal distribution of brain neural flexibility was also revealed. Finally, we defined the brain regions with significantly higher neural flexibility than that of whole brain as the flexible club and compared them to the diverse club (44) and functional hub (4).

Diverse club was defined as brain regions exhibited a high participation coefficient (PC). The PC indicates how evenly a node’s connections are distributed to different modules. For a given modular assignment, the participation coefficient of node can be calculated as

where is the sum of node ’s edge weights, is the sum of node ’s edge weights to module , and is the total number of modules. Since PC calculation requires modular assignments, we performed the (classic) Louvain algorithm (71) for each subject.

Brain functional hub was defined as the brain regions exhibiting high BC. Given a graph, BC is defined as the fraction of shortest paths between node pairs that travels through the node of interest:

where is the total number of shortest paths between node and node within the graph , is the total number of these shortest paths passing through node , and is the total number of nodes in the graph.

To determine the unique features associated with the flexible club when compared to the better-documented functional hub and diverse club, spatial overlapping ratio, mean/SD of dynamic connectivity strength, and resilience to target attack were compared. Specifically, for each region and time window, the mean connectivity strength to all other regions was calculated. The mean (i.e., averaging across all time windows) and SD of dynamic connectivity strength were calculated for each region and further spatially averaged across regions for each of the three club types (i.e., flexible club, diverse club, and functional hub). In addition, brain targeted attack was simulated by removing regions (i.e., flexible club, diverse club, and functional hub) and all their connections, and brain global efficiency (4) was calculated. The global efficiency of node can be calculated as

where is the number of nodes of the whole graph and is the minimum path length between the node and all other nodes in the whole graph. The mean GE was obtained by averaging across all nodes. A paired t test was performed to compare the statistical difference between these parameters.

In this study, statistical significance was considered as . For multiple comparison correction, FDR correction was performed by controlling FDR at a level .