Spatiotemporal cerebral blood flow dynamics underlies emergence of the limbic-sensorimotor-association cortical gradient in human infancy

Ouyang, Minhui; Detre, John A.; Hyland, Jessica L.; Sindabizera, Kay L.; Kuschner, Emily S.; Edgar, J. Christopher; Peng, Yun; Huang, Hao

doi:10.1038/s41467-024-53354-7

Download PDF

Article
Open access
Published: 17 October 2024

Spatiotemporal cerebral blood flow dynamics underlies emergence of the limbic-sensorimotor-association cortical gradient in human infancy

Minhui Ouyang ORCID: orcid.org/0000-0001-8013-2553^1,2,
John A. Detre^2,3,
Jessica L. Hyland¹,
Kay L. Sindabizera¹,
Emily S. Kuschner^1,4,
J. Christopher Edgar^1,2,
Yun Peng⁵ &
…
Hao Huang ORCID: orcid.org/0000-0002-9103-4382^1,2

Nature Communications volume 15, Article number: 8944 (2024) Cite this article

3510 Accesses
4 Citations
1 Altmetric
Metrics details

Subjects

Abstract

Infant cerebral blood flow (CBF) delivers nutrients and oxygen to fulfill brain energy consumption requirements for the fastest period of postnatal brain development across the lifespan. However, organizing principle of whole-brain CBF dynamics during infancy remains obscure. Leveraging a unique cohort of 100+ infants with high-resolution arterial spin labeled MRI, we find the emergence of the cortical hierarchy revealed by the highest-resolution infant CBF maps available to date. Infant CBF across cortical regions increases in a biphasic pattern featured by initial rapid and subsequently slower rate, and break-point ages increasing along the limbic-sensorimotor-association cortical gradient. Increases in CBF in sensorimotor cortices are associated with enhanced language and motor skills, and frontoparietal association cortices with cognitive skills. The study discovers emergence of the hierarchical limbic-sensorimotor-association cortical gradient in infancy and offers standardized reference of infant brain CBF and insight into the physiological basis of cortical specialization and real-world infant developmental functioning.

Assessing cerebral blood flow, oxygenation and cytochrome c oxidase stability in preterm infants during the first 3 days after birth

Article Open access 07 January 2022

Diffuse correlation spectroscopy blood flow monitoring for intraventricular hemorrhage vulnerability in extremely low gestational age newborns

Article Open access 27 July 2022

Cognitive and psychiatric relevance of dynamic functional connectivity states in a large (N > 10,000) children population

Article Open access 31 July 2024

Introduction

Although human brain growth continues into adulthood, it is most dramatic in infancy, doubling in size during the first year and to almost adult size by age 5^1,2. The brain is one of the most metabolically expensive organs in the body; despite representing only about 2% of total body weight, the adult human brain consumes nearly 20% of the body’s basal metabolism³. However, the infant brain uses a strikingly larger fraction of the total body energy consumption, exceeding 40%^4,5. Such high energy consumption is vital for supporting intensive maturational processes, including synaptogenesis, glial proliferation, and myelination⁶. Cerebral blood flow (CBF) delivers both oxygen and glucose, the principal brain substrates for energy production, to the brain. Because of the tight coupling between regional CBF (rCBF) and local cerebral metabolism^7,8, rCBF can be used as a surrogate measure of regional brain metabolism. The “brain growth spurt” during infancy⁹ is associated with enhanced vulnerability to disorders¹⁰. If blood flow is interrupted or inadequately supplied, the brain is more susceptible to disorders such as ischemic stroke and hypoxic ischemic brain injury. Infancy is also a crucial period of behavioral development, with proper rCBF needed to support the emergence of cognitive abilities and behaviors. Accordingly, characterization of normative infant rCBF trajectories provides an opportunity to explore how brain regions use energy differentially in support of their functional specializations in health and disease, as well as the physiological basis for the progression of real-world developmental functioning. However, information on how cortical rCBF topography evolves throughout infancy and its associations with behavioral development is surprisingly lacking.

Cortical hierarchy, representing a more pronounced form of the differential cortical development along specific gradients^11,12, is an organizing principle of human brain development. Although the sensorimotor-association gradient has been well characterized in adolescents and adults¹³, such a pattern cannot be completely discerned at birth. During infancy, cortical gradient may exhibit a slightly different pattern than in later developmental stages. The gradient could involve other components such as the limbic system, given the phylogenetically old nature of limbic functions. Tremendous developmental brain changes occurring throughout infancy are critical for the emergence of cortical gradients. Neuroimaging derived measures of infant brain¹⁴ such as cortical morphology e.g., refs. ^15,16, white matter microstructure e.g., refs. ^2,17,18,19^., and functional connectivity e.g., refs. ^20,21,22 suggest a hierarchical sequence of cortical development in which limbic and primary areas mature before association cortices. These developmental brain changes demand rapid increases in metabolic substrate. Previous fluorodeoxyglucose positron emission tomography (FDG-PET) studies suggest a spatially heterogenous courses of metabolic maturation across infant cortex, with increases in local cerebral glucose utilization observed in temporal and primary visual cortices at 2–3 months of life, and then in the frontal cortex at 6–12 months^23,24. The dynamic course of local brain metabolism revealed in these studies provides grounds for anticipating a spatiotemporally heterogenous pattern of rCBF trajectory during infancy. We thus hypothesized that the infant rCBF trajectory proceeds in a hierarchical manner, beginning with limbic and primary cortices at the bottom and gradually transitioning to heteromodal association areas toward the top of the hierarchical levels.

To date, there has been little, if any, direct rCBF-based physiological evidence supporting the emergence of limbic-sensorimotor-association cortical gradient. Complicated by both neural and non-neural components^25,26, the mechanisms of fMRI blood-oxygen-level-dependent (BOLD) signal remain elusive, and resultantly fMRI-based mapping of cortical hierarchy likely reflect measures of neural activities mixed with other factors such as vasomotor activity. Complementary to the valuable fMRI BOLD studies, our study exclusively based on fundamental physiological CBF measurements that directly and quantitively reflect regional resting state neural activity, could offer refreshing insight into the emergence of the limbic-sensorimotor-association cortical gradient. Arterial spin labeled (ASL) MRI uses magnetically labeled arterial blood water protons as an endogenous tracer for reliably estimating rCBF²⁷. As compared to FDG-PET, which requires the intravenous administration of a radioactive tracer, ASL MRI is completely noninvasive and hence well suited for use in healthy infants and children e.g., refs. ^{28,29,30,31,32}. Collectively, there is still a knowledge gap in not only mapping infant rCBF to regional precision but also depicting the organizing principle of rCBF dynamics that underpin cortical gradient emergence during the critical developmental period from 0 to 2 years.

In this study, by acquiring the highest-resolution infant CBF maps available to date with densely sampled ages during infancy, we explored the organizing principle of spatiotemporal dynamics of rCBF in typically developmental infants with greater precision. The high-resolution spatiotemporal dynamics of rCBF were obtained by capitalizing on a state-of-the-art pseudo-continuous ASL (pCASL) imaging protocol with 3D multi-shot stack-of-spirals readout^33,34 and background suppression. Unlike previous studies e.g., refs. ^30,31 in the setting of clinical suspicion of diseases, a typically developmental infant CBF study with such a scale of samples within 0–2 years offers invaluable early typically developmental physiological dataset. We found that the global CBF increased in a logarithmic fashion during infancy. The high-resolution, high-quality infant rCBF maps revealed physiological heterogeneity in spatial distribution as well as age-related increases. These findings were also replicated in a second cohort using high-resolution 2D pCASL imaging. By employing a data-driven clustering approach to identify topological regionalization of rCBF with a similar changing pattern, we further discovered three distinct clusters, including limbic, sensorimotor, and frontoparietal areas that reflect the limbic-sensorimotor-association cortical gradient. Lastly, we found that rCBF changes were associated with improved developmental functioning across multiple domains. Together, these findings suggest that the spatiotemporal infant rCBF dynamics underlie the emergence of a hierarchical limbic-sensorimotor-association gradient and contribute to the development of behavior in infancy, setting the stage for cortical function specialization.

Results

Infant participant profile

A total of 134 infants from two cohorts (N = 78 for cohort-1 and N = 56 for cohort-2; see “Participants” in “Methods” section and Supplementary Table 1) participated in this study. Briefly, cohort-1 is from an ongoing study at the Children’s Hospital of Philadelphia. This study investigates normal brain development during infancy. All infants from cohort-1 are typically developing children and recruited completely for research purposes. Cohort-1 infant MR images were collected during natural sleep. Cohort-2 is from Beijing Children’s Hospital, where cohort-2 infants were clinically referred for neurological pathology, and cohort-2 infant MR images were collected under sedation (see “Participants” in “Methods” section for details). Each infant in cohort-2 had their clinical history thoroughly reviewed to rule out developmental abnormalities, and no brain abnormalities were detected through reading of these infants’ MRI scans by an experienced pediatric radiologist (YP). It should be noted that except that Fig. 1 was from both cohorts, almost all results (Figs. 2–6 and almost all Supplementary Figs.) were obtained from one cohort (cohort-1) of infants in one MR acquisition ___location.

**Fig. 1: Developmental curve of global cerebral blood flow (CBF) throughout infancy.**

**Fig. 2: Precise physiological variability at finer detail across infant brain regions.**

**Fig. 3: Nonuniform age-related increases of regional CBF during infancy.**

**Fig. 4: Infant rCBF increases according to a hierarchical limbic-sensorimotor-association gradient.**

**Fig. 5: Regionally specific rCBF increases are associated with developmental functioning in infants.**

**Fig. 6: Infant rCBF reorganized by spatiotemporally varying energy demand for brain maturation.**

Logarithmic increase of global CBF during human infancy revealed by datasets of 100+ infants with densely sampled ages

As a first step, we established the infant developmental curve of global CBF with a relatively large datasets from 119 infants (60M/59F; Age range of 1.4–28 months with mean/standard deviation 12.5 ± 7.1 months) who had completed phase-contrast (PC) perfusion MRI in both cohort-1 and cohort-2 (see “Participants” in “Methods” section and Supplementary Table 1). Global CBF for each infant was measured by PC perfusion MRI, a reliable global CBF measurement technique³⁵ that encodes the velocity of flowing spins as phase changes in an MRI image. By measuring the flow flux at the four feeding arteries of the brain (Supplementary Fig. 1a), including bilateral internal carotid arteries and vertebral arteries, global CBF can be noninvasively measured. We investigated the effects of age, sex, and cohort on global CBF using a large infant database from these 119 infants (see “Participants” in “Methods” section). We found that global CBF increased in logarithmic manner during infancy (r = 0.823, P < 2 × 10⁻¹⁶; Fig. 1), after testing and quantifying the relative goodness-of-fit among different models (including linear, logarithmic, exponential, Poisson, and quadratic polynomial) for fitting the developmental curve (Supplementary Fig. 2a and Supplementary Table 2). Sex had no effect on global CBF trajectory, nor did cohort have any detectable effect (Supplementary Fig. 2b and Supplementary Table 2). Segmented regression analysis revealed two distinct phases of global CBF increases during infancy (F (2,117) = 22.23, P < 6.88 × 10⁻⁹; Supplementary Table 4), with a similar breakpoint as that of the total daily energy expenditure growth³⁶ (see “Discussion” section). Global CBF rose rapidly at 3.52 ± 0.45 ml/100 g/min per month in the first phase from birth to a break point at 10.75 months of age (95% CI—confidence interval: 8.70, 12.81), followed by sequentially slower increases at a rate of 0.43 ± 0.23 ml/100 g/min per month in the second phase. These findings highlight significant physiological changes that occur during infant brain development and offer a normative reference of global CBF during infancy.

Regional physiological variability at finer detail across infant brain revealed by high-resolution rCBF maps based on arterial spin labeled (ASL) perfusion MRI

Having established the normative developmental trajectory of global CBF, we then examined how regional CBF topography evolves with age during infant development. Due to the small size of infant brains, substantially higher spatial resolution usually at the cost of signal-to-noise ratio (SNR)² of rCBF images is required to achieve comparable “anatomic” details to the adult brain. To achieve high spatial resolution and high SNR, we leveraged an advanced 3D multi-shot, stack-of-spiral pCASL MRI sequence^33,34 and developed an infant-dedicated protocol to acquire in vivo ASL images at 2.5 mm isotropic resolution. In contrast to resolutions as low as 4–10 mm conventionally used in infant and pediatric populations, the resulting datasets acquired at the single Philadelphia site (N = 76, 32M/44F; Supplementary Table 3) provided the highest resolution ASL data to date for infants, and revealed heterogeneous spatial distribution patterns (Fig. 2 and Supplementary Fig. 3) that could not be observed previously. Higher resolution and SNR with the cutting-edge pCASL^33,34 can be clearly observed when compared with conventional protocol (Supplementary Fig. 4). The dense acquisition of thin axial slices in our infant-dedicated pCASL protocol facilitates more precise estimation of regional CBF. These highest resolution infant CBF maps available to date capture around 3 times of cortical voxels captured in the conventional rCBF maps (Supplementary Fig. 4), revealing the unprecedented details of the whole-brain physiological dynamic patterns in infants. The rCBF findings shown in Figs. 2–6 are from 76 typically developing infants in cohort-1 only. To ensure reliability of the datasets, we examined the reproducibility of rCBF maps (Supplementary Fig. 1b) by scanning an infant twice using the same scanner and identical pCASL protocol. Consistent with literature on high test-retest reliability of ASL MRI³⁷, the intraclass correlation coefficient of two scans in this study was 0.869 (95% CI: 0.861,0.876) (Supplementary Fig. 1c), indicating excellent reliability of rCBF maps obtained with our pCASL protocol.

The resulting higher resolution and effective SNR of infant rCBF maps enable more precise localization and provides finer details for physiological variability in the infant brains than prior work^29,30,31. As demonstrated by rCBF maps from three representative infants aged 2, 8 and 17 months in Supplementary Fig. 3, the densely acquired axial slices of the infant rCBF maps reveal a high gray and white matter contrast (Supplementary Fig. 3a–c). After mapping to the cortical surface, regional heterogeneity of rCBF across infant cortex can be clearly appreciated (Supplementary Fig. 3d–f). Specifically, we found that rCBF values in the 2-month-old infant brain are higher in primary cortex (sensorimotor, auditory, and visual) and lower in association cortex (frontal and parietal), whereas in the 17-month-old infant, higher rCBF values are more prominent in association cortex. We further calculated group-averaged rCBF maps at six milestone age groups (Fig. 2). Consistent with the spatial distributions observed in individual infants, higher rCBF values in group-averaged rCBF maps shown in Fig. 2 can be found in primary sensorimotor, auditory and visual cortices, indicated by white, purple and orange arrows, respectively, in 0–3-month infants versus association cortex in infants older than 12 months. These observations demonstrate good spatial resolution and SNR of the infant rCBF maps to capture spatiotemporal variability for us to further tackle the emergence of cortical hierarchy from the physiological aspect.

Nonuniform age-related regional CBF increases across heteromodal and unimodal cortex during infancy

Based on the empirical observation that higher rCBF values are located in the unimodal primary cortex of younger infants and in the heteromodal association cortex of older infants (Fig. 2 and Supplementary Fig. 3), we hypothesized that age-related changes in rCBF would vary heterogeneously across cortical regions. To test this hypothesis, we used the high-resolution infant rCBF dataset from cohort-1 (N = 76) and logarithmic models to examine how rCBF is associated with age during infant development, based on a logarithmic pattern for global CBF shown in Fig. 1. Sex and in-scanner motion were included as covariates in these models. After correcting for multiple comparisons with the Bonferroni method (Z ≥ 5.1, Bonferroni-corrected p value P_bonf ≤ 0.05), we found that age-related increases in rCBF were nonuniform across the cortex (Fig. 3a). Specifically, the most significant age-related increases in rCBF were seen in areas of heteromodal association cortices, including prefrontal cortex, inferior parietal lobule, lateral temporal cortex, precuneus, and posterior cingulate cortex. Substantial rCBF increases were also observed in unimodal sensorimotor cortex, but to a lesser extent.

Our findings were further supported by a validation analysis of voxel-wise age effects using rCBF maps (N = 48) obtained from cohort-2 (Beijing site)³⁰. The age effects were highly reproducible in both infant rCBF datasets, showing similar spatial topographies (Supplementary Fig. 5a, b) and high between-dataset correspondence (r = 0.426, P_perm = 0.0001; Pearson correlation with permutation test; Supplementary Fig. 5c). The infant rCBF dataset from cohort-1 in Philadelphia (Fig. 3a and Supplementary Fig. 5a) tended to yield stronger age effects than the validation cohort from Beijing (cohort-2; Supplementary Fig. 5b), most likely due to differences in pCASL acquisition (voxel size of 2.74 × 2.74 × 5 mm³ for the Beijing cohort versus the higher resolution of 2.5 × 2.5 × 2.5 mm³ for the Philadelphia cohort). Additionally, infants in the Beijing cohort (cohort-2) were scanned due to clinical indications, although no brain abnormalities were detected in cohort-2 infants based on an experienced pediatric radiologist’s MRI reading. Infants in the Philadelphia cohort (cohort-1) were typically developmental ones recruited solely for research.

We further examined the age effects of regional CBF using both predefined functional networks from Yeo brain network atlas³⁸ (Fig. 3b, c) and a data-driven clustering approach (Fig. 4) with high-resolution rCBF dataset from the cohort-1. Notably, age effects were nonuniformly distributed across seven functional networks (Fig. 3b), with the largest increases in rCBF occurring in fronto-parietal and default-mode networks, and with modest increases in rCBF observed in sensorimotor and limbic networks. For instance, rCBF from the fronto-parietal and default-mode network regions demonstrated stronger age effects (FPN-rep: β = 24.08, r = 0.891, P < 2 × 10⁻¹⁶; DMN-rep: β = 20.21, r = 0.887, P < 2 × 10⁻¹⁶) than rCBF from sensorimotor and limbic network regions (SM-rep: β = 15.44, r = 0.757, P = 6.51 × 10⁻¹⁴; LIM-rep: β = 16.59, r = 0.632, P = 6 × 10⁻⁹) (Fig. 3c). We constructed rCBF maps from the developmental curves at various ages throughout infancy (Supplementary Figs. 6–7). The resulting trajectory was animated to create a quantitative time-lapse movie of infant rCBF dynamics (Supplementary Movie 1). Together, the spatiotemporal rCBF measures revealed a shifting pattern of distribution with its developmental gradient beginning in limbic and sensorimotor areas and spreading rostrally into prefrontal cortex, as well as dorsally and caudally into parietal and temporal cortices.

Hierarchical limbic-sensorimotor-association gradient representing a general organizing principle of infant regional CBF dynamics

The findings above demonstrate notable areal heterogeneity in the age effect of rCBF, prompting the brain to be divided into clusters (e.g., groups of voxels), each with a more similar rCBF increase pattern than the rest of the brain. To provide a more nuanced understanding of the physiological dynamics during infancy, we performed data-driven clustering analysis of infant rCBF maps using non-negative matrix factorization (NMF)³⁹ to identify distinct clusters. The NMF approach has been shown to identify brain networks robustly and accurately using neuroimaging data e.g., ref. ¹¹. Only high-resolution rCBF maps from 76 TD infants from cohort-1 were included in these data-driven clustering analyses. We discovered three clusters along the hierarchical limbic-sensorimotor-association gradient with unique rCBF increase patterns (Fig. 4a and Supplementary Fig. 8). By comparing with the commonly used Yeo brain network atlas³⁸, voxels corresponding to each cluster primarily highlighted limbic, sensorimotor, and frontoparietal regions, reflecting limbic/sensorimotor-association hierarchy. To further examine the break-point ages between faster and slower rCBF increases, we conducted segmented regression analysis (elaborated in global CBF analysis) on averaged rCBF values from each cluster, resulting in similar biphasic patterns for all clusters (all F (2,74) ≥ 5.61, p ≤ 5.4 x 10⁻³). RCBF measures for all clusters increase faster in the first phase and then increase steadily but at a much slower rate in the second phase (Fig. 4a and Supplementary Table 4). The identified break points differed across clusters and exhibited a clear pattern along the hierarchical gradient. Notably, the break point was earliest for the limbic cluster at 6.71 months (95% CI: 4.09, 9.33), followed by the sensorimotor cluster at 7.69 months (95% CI: 5.63, 9.75), and then lastly the frontoparietal cluster at 10.36 months (95% CI: 8.21, 12.50). We also examined the break point profiles of all voxel-wise rCBF trajectories in the infant cortex. The resulting histograms (Fig. 4b) demonstrate that voxels in cortex from higher hierarchical levels (e.g., frontoparietal regions) have older break-point ages than voxels from lower hierarchical levels (e.g., sensorimotor and limbic regions). The robustness of our findings at the cluster and voxel-wise levels demonstrates that the biphasic pattern of rCBF increases varies along a cortical hierarchy, from limbic, and sensorimotor to frontoparietal association cortex.

Rates of regional infant rCBF increase calculated as the first derivative of fitted logarithmic curves of all cortical voxels further supports the limbic-sensorimotor-association gradient. We observed the highest rates of rCBF increase before 6 months, followed by dramatically slower increase rates after 6 months (Fig. 4c). Intriguingly, the duration of the rapid increase is longer in regions with higher hierarchical levels across three clusters (Fig. 4c). The rCBF increase rate maps at six milestone ages are demonstrated in Fig. 4d.

Regionally specific rCBF increases supporting the emergence of infant behavior and developmental functioning

With the NMF parcellation in Fig. 4a, we sought to determine where the rCBF increases support certain enhanced behavior during infancy by investigating the association between the rCBF and measures of infant developmental functioning. Using the Bayley scales of infant and toddler development, we assessed infants’ behavior and functioning in the areas of motor, language, and cognition. We employed generalized additive models (GAMs)⁴⁰ to model the association between all voxel-wise rCBF measures and infant neurodevelopmental outcomes, controlling for age, sex, in-scanner motion, and family socioeconomic status. Among the 76 infants who had high-resolution pCASL acquired in cohort-1, 49 of them have their Bayle scales and family socioeconomic status available and are included in the association analysis here. Cohort-2 is not included for the rCBF-behavior analyses. The results were clusters of voxels significantly correlated with Bayley ___domain scores (Fig. 5a–c). Specifically, the largest clusters of infant rCBF significantly correlated with motor (e.g., R.ANG, t = 3.670, P_bonf = 0.007; Fig. 5a) and language skills (e.g., L.MOG/IOG, t = 3.703, P_bonf = 0.006; Fig. 5b) were found in the sensorimotor/frontoparietal cortices. RCBF cluster located in the frontoparietal cortex (e.g., R.ANG, t = 3.624, P_bonf = 0.007; Fig. 5c) significantly correlated with cognitive skill. The complete list of clusters is shown in Supplementary Table 5.

We further examined how behavioral manifestations were encoded by the identified hierarchic limbic-sensorimotor-association gradient. As shown in Fig. 5d, the infant rCBF hierarchic distribution regions were largely associated with a specific developmental functioning. In particular, language skill was primarily encoded in sensorimotor cortices, and cognitive was mainly encoded in frontoparietal cortex. Motor skill was encoded in sensorimotor and frontoparietal cortices. These findings suggest that region-specific rCBF increases along the hierarchic limbic-sensorimotor-association gradient are associated with enhanced developmental functioning during infancy.

Infant rCBF reorganized by spatiotemporally varying energy demand for brain maturation

Given that blood flow delivers glucose and oxygen to different brain regions to meet energy demands, rCBF distributions are tightly associated with the neuronal activities and functional specializations among brain areas⁷. Consistent with the literature^7,8, and shown in Fig. 6a, rCBF is a high-fidelity surrogate measure of cerebral metabolic rate of glucose (CMRglc) measured by FDG-PET. Despite individual differences, a strong spatial association between adult rCBF and CMRglc maps (r = 0.594, P_perm = 0.0001; Pearson correlation with permutation test) was observed (Fig. 6a). However, FDG-PET is usually not applicable to typically developing infants given radiation concerns. From Fig. 6b, infant rCBF is progressively reorganized during development from significantly different distribution at a younger age (e.g., 1.5 months) indicated by lower r value (r = −0.035, P_perm = 0.0001) to a distribution highly aligned with the adult metabolic distribution map at an older age (e.g., 28 months) indicated by higher r value (r = 0.635, P_perm = 0.0001). A reorganization plateau is reached at around 10 months (Fig. 6b), coinciding with the break-point age of global CBF increase (Fig. 1). These findings suggest that the cerebral cortex undergoes a dynamic course of metabolic maturation during infancy, and the hierarchy-dependent age-related changes of rCBF are impacted by the allocation of brain energy demand for brain maturation.

Discussion

We found the emergence of the limbic-sensorimotor-association cortical hierarchy in the infant brain using noninvasive physiological imaging based on cerebral blood flow. This physiological cortical hierarchy reflects changes in regional neural development that likely underlies cortical functional specialization. We also discovered the physiological basis of behavioral manifestation by associating increases in cerebral blood flow in specific cortical regions and corresponding behavioral measures. Our established normative charts of rCBF at the highest resolution available during infancy not only advances current understanding of normal brain development but also provides an anchor point for detecting abnormal brain perfusion (e.g., brain injury) during this vulnerable period. Taken together, we filled the knowledge gap by elucidating the rCBF dynamics underpinning emergence of the critical hierarchical limbic-sensorimotor-association cortical gradient and revealing rCBF-behavior association in infancy through acquiring the highest-resolution infant rCBF available to date.

Maintaining proper CBF dynamics has a profound impact on brain development and function. CBF delivers glucose and oxygen, the main substrates for energy metabolism in the brain. It is closely coupled with the cerebral metabolic rates for glucose (CMRGlu) and oxygen (CMRO2)^7,8, the direct measures of the rate of energy consumption. Loss of CBF is accompanied by irreversible damage to neurons and synapses⁴¹. The massive and rapid development renders infancy vulnerable to disorders^9,10. Thus, infant brain is particularly susceptible to cerebral ischemia and stroke, which disrupts CBF, leaving the brain with inadequate energy to support its underlying maturational processes. The developmental pattern of global CBF during infancy increases rapidly and logarithmically with age and coincides with known dynamic changes in brain metabolism^23,24, reflecting a phase of escalating metabolic demands associated with the overproliferation of energetically costly synapses and dendritic spines⁴². Dramatic cellular processes in early developmental brains, such as extensive synaptogenesis, dendritic arborization, glial proliferation, and myelination, are critical for the formation and growth of neural circuits⁶, allowing for rapid brain function development. These intensive developmental processes and increased neural activity necessitate higher metabolic demands, promoting angiogenesis for an elevation in CBF to supply more glucose and oxygen. Moreover, we found that global CBF had a similar rise and break point as the adjusted total daily energy expenditure during infancy³⁶, consistent with the notion that much of infant whole-body metabolism is dedicated to brain development⁵.

We further observed that rCBF was heterogeneously distributed across brain regions during infancy. The spatial pattern of infant rCBF resembled that of local brain metabolism at a certain age. For example, the prominent physiological activity in primary sensorimotor, auditory and visual cortical areas, putamen, and thalamus in 0-3-month-old infants is consistent with higher glucose utilization^23,24 in these regions than the rest of the brain. In contrast, by 9 to 12 months, infants’ rCBF patterns qualitatively appear like those of adults, with strong physiological activity in frontal and parietal association cortices, matching the changes in their metabolic patterns as well. We also discovered that infant rCBF increased in a spatially heterogenous fashion, resulting in a shift of rCBF distribution patterns as a function of age. These findings align with prior infant ASL research e.g., refs. ^29,30,31, which noted uneven distributed rCBF and its changing patterns but mainly focused on limited brain regions and lacked sufficient coverage of the entire first 2 years of life. The spatiotemporal heterogeneity of infant rCBF dynamics underlies the emergence of cortical gradient.

Most notably, our study provides the first clear view of physiological hierarchy along the limbic-sensorimotor-association gradient throughout infancy. To our knowledge, such cortical axis with limbic component has not been reported so far, although the sensorimotor-association gradient been proposed for cortical maturation at much later developmental stages such as adolescence and adulthood e.g., ref. ⁹. The identified physiological hierarchy aligns well with the pronounced and qualitatively distinct regional differences in infant functional connectivity^20,21,22,30, neuroanatomy^2,16,43, and transcriptomic profile^44,45. For instance, limbic and primary sensorimotor functional networks are relatively well established at birth^20,21. In contrast, functional networks involving association cortices (e.g., default mode, frontoparietal, and attention networks) are often recognizable in topographically incomplete forms at birth and then grow rapidly during the first two years of life^22,30, gradually establishing the limbic-sensorimotor-association functional gradient. Thus, it is not surprising that rCBF rises greatly in these association cortices to meet the extra energy demand for functional network and gradient formation³⁰. Infant brain also undergoes remarkable but non-uniform morphological changes, with lateral temporal, parietal, and frontal cortices expanding more quickly than primary sensorimotor and occipital cortices¹⁶. Cortical thickness of the medial orbitofrontal cortex from limbic network reaches its peak earliest with respect to all cortical cortex¹⁶. White matter develops in a similar spatial pattern^2,12. For example, limbic fibers (e.g., fornix) appear and myelinate earliest, followed by sensorimotor projection fibers (e.g., corticospinal tract), and then association fibers connecting to frontal areas (e.g., inferior fronto-occipital fasciculus)^2,17,18,19. These rapid functional and structural brain changes are reliant on the spatiotemporal regulation of the transcriptome with gene expression during infancy showing differences between primary and association cortical areas^44,45 and pronounced enrichments of biological processes related to synaptogenesis and myelination that proceed heterochronically across brain regions⁴². The emergence of the limbic-sensorimotor-association gradient comes at substantial metabolic cost. Developmental refinement of rCBF dynamics throughout infancy reflects the reallocation of the brain’s energy resources to promote the establishment of this cortical gradient. This infant cortical hierarchy with limbic component may transition over time, merging into the sensorimotor-association gradient observed in adolescents and adulthood. Our findings converge with the functional, transcriptomic and structural findings in the literature and fill a gap regarding the physiological hierarchical emergence of the human cortex. The CBF-based physiological hierarchy may underlie the emergence of macrostructural, microstructural, and functional hierarchy.

In parallel with the rapid brain changes for establishing cortical hierarchy, infants begin to walk, talk, and develop increasingly complex behavior. We found that increased rCBF was associated with measures of infant behaviors and developmental functioning. As the majority of brain’s energy is devoted to neural signaling processes that maintain brain function⁴⁶, increased CBF may ensure that the brain meets the energy demands needed for the formation of large-scale functional networks capable of supporting complex behaviors^30,47. Cortical regions with significant rCBF-behavior associations exhibited selectivity for different behaviors, suggesting distinct brain circuits facilitating the development of specific behaviors and functioning. For language development, increased rCBF was predominantly seen in sensorimotor regions that contributing to the execution of motor activities. This is in line with the finding of sensorimotor network influences speech and language acquisition in infancy⁴⁸. In contrast, cognitive development was more associated with rCBF in frontoparietal areas. It is noteworthy that the frontoparietal network plays an important role in infants’ social cognitive development⁴⁹. The observed association of several cortical regions with multiple Bayley subscales likely reflects the inherent intercorrelations⁵⁰ within the subscales due to similar procedures utilized across domains or overlap in skills required to complete subtests. Understanding how rCBF dynamics support the emergence of infant behavior and developmental functioning may provide potential biomarkers for early detection of disorders in infants, as adaptive behavior delays during infancy have been linked to later diagnoses of neurodevelopmental conditions such as autism⁵¹.

Prior studies that focused on small subsets of brain regions or incomplete coverage of the first 2 postnatal years found varied rCBF distribution and changes in infant brains e.g., refs. ^29,30,31, but our study is the first to delineate the rCBF hierarchy throughout infancy. Low-resolution maps in prior work hampered the capacity to precisely chart the infant rCBF topography. Capitalizing on a cutting-edge 3D multi-shot, stack-of-spiral pCASL with background suppression^33,34, the much higher imaging resolution in present study distinguishes rCBF heterogeneity across proximate brain areas. The resultant highest infant rCBF maps to date with 2.4 times or more resolution improvement than previous studies change the paradigm and enable capturing the unprecedented details of the infant whole-brain physiological hierarchy throughout infancy. According to the ASL consensus paper²⁷, 3D acquisition is the recommended ASL readout method because it is SNR-efficient and can be optimally combined with background suppression to effectively reduce physiological noise when compared to 2D readout. In this study, only 3D ASL data from cohort-1 were used for the main results, whereas 2D ASL data from cohort-2 was used for validation only. Due to varying brain sizes and cortical thicknesses across infants, using a consistent pCASL acquisition resolution could lead to heterogeneous partial volume effects. The biased measures of rCBF due to partial volume effects were alleviated by high-resolution 3D ASL acquisition and were further minimized with individual PC MRI calibration³⁰. Partial volume effects can also be alleviated by measuring rCBF values at the core “skeleton” e.g., ref.¹¹ of the cortical ribbon given robust infant cortical segmentation in the future studies. Furthermore, we mitigated the impact of in-scanner head motion on our results by image scrubbing to remove image volumes with large motion, rigorously assessing image quality for data inclusion and co-varying for small motion in all analyses. These procedures limit the likelihood that the reported findings are due to in-scanner head motion. Despite slightly more female in cohort-1, sex was systematically accounted for as a covariate in all analyses to minimize any potential effects. We demonstrated that rCBF is strongly correlated with local brain metabolism at rest despite that rCBF and CMRglc were acquired from different adult individuals. The observed age-related rCBF change pattern in infants of this study is also aligned with CMRglc measures from previous studies^23,24 However, caution needs to be taken for the interpretation as rCBF is not a direct measure of the rate of brain energy consumption. With only cross-sectional data presented in this study, longitudinal studies are warranted to minimize the inter-subject variability when delineating age-related rCBF dynamics in infancy. Future studies can investigate how combining rCBF with neuroimaging measures from other modalities, such as functional MRI and diffusion MRI, might reveal multidimensional principles of infant brain maturation. This work is a step towards the normative rCBF milestone for clinical translation. With rapid MR technical advances in hardware and sequence, the high-resolution rCBF currently obtained with a research protocol can be soon translated into clinical settings. It is recognized that a barrier to ASL clinical translation is the lack of standardization of these perfusion MRI methods across scanners and across sites. With the growing effort on generating consensus e.g., refs. ^27,52 and promoting standardization for imaging biomarker (e.g., QIBA: Quantitative Imaging Biomarkers Alliance from the Radiological Society of North America, https://qibawiki.rsna.org), variation of rCBF measurements across vendors and sites will likely be reduced over time as more studies incorporate ASL and more clinical data are collected and validated. Besides other imaging modalities such as functional MRI, perfusion MRI adds a new avenue of tackling behavioral manifestations from an imaging point of view. Future studies with larger cohorts and potentially more specific behavioral assessments are warranted to comprehensively characterize the rCBF-behavior relationship in infants.

Methods

Participants

As shown in Supplementary Table 1, a total of 134 infants from two independent cohorts were recruited in this study. In cohort-1 (N = 78), we collected high-resolution pCASL of isotropic 2.5 mm and PC perfusion MRI in infants at the Children’s Hospital of Philadelphia. In cohort-2 (N = 56), infants underwent pCASL and PC perfusion MRI at the Beijing Children’s Hospital. We combined PC-MRI data from two cohorts to investigate the developmental model of global CBF in a large sample (Fig. 1; N = 119; 63 from cohort-1 and 56 from cohort-2). No cohort effects were observed (Supplementary Table 2). We used high-resolution pCASL data from cohort-1 (N = 76) to delineate the developmental models of regional CBF during infancy (Figs. 2–6). For validation, key findings reported for voxel-wise age effects on regional CBF during infancy were replicated using the pCASL data collected from infant cohort-2 (N = 48; Supplementary Fig. 5). Below are detailed descriptions of each infant cohort.

Cohort-1 is part of an ongoing study investigating typical brain development during infancy at the Children’s Hospital of Philadelphia, Philadelphia, USA. Participants were recruited from the greater Philadelphia metropolitan area with assistance from the Clinical Reporting Unit at Children’s Hospital of Philadelphia for identifying eligible participants and through study flyers. All recruited infants were enrolled as typically developmental infants (i.e., brain imaging not clinically indicated). They were recruited completely for research purposes, which involved studying normal infant brain development. All infants were selected after rigorous screening procedures conducted by an experienced psychologist (ESK) based on their medical history and family medical history. Exclusionary criteria for infants were: premature birth or significant birth complications; known neurological disorders; genetic conditions with known neurological or neurodevelopmental sequelae; severe tics, or severe head trauma; known developmental or medical diagnoses; major congenital heart disease or congenital infection; uncorrectable sensory (hearing, visual) issues; having a sibling or parent with a diagnosis of autism spectrum disorder or schizophrenia; and 1.5 or more standard deviations below the mean on the behavioral assessment of developmental milestones. From the 275 infants who were originally recruited to a larger normal brain development study, 78 infants who completed scans of either PC-MRI or pCASL perfusion MRI participated in this study as cohort-1. Out of 78 infants, advanced pCASL was successfully acquired from 76 infants aged 1.2–28 months (mean/standard deviation 10.98 ± 6.82 months; 32M/44F), and PC-MRI was acquired from 63 infants aged 1.5–28 months (mean/standard deviation 9.84 ± 7.36 months; 26M/37F). MRI scans of infants in cohort-1 were collected during their natural sleep. Perfusion MRI sequences (pCASL and PC-MRI) were among the last in the almost 1-hour long MRI protocol. This led to a relatively small number of completed pCASL datasets compared to the total number of infants recruited, as many infants woke up before the pCASL could be acquired. All procedures for cohort-1 were approved by the Institutional Review Board at the Children’s Hospital of Philadelphia in compliance with ethical regulations and standards (Approval number IRB16-013203). Written informed consent was obtained from all infants’ parents.

Cohort-2 included 56 infants aged 1.4–27.7 months (mean/standard deviation 14.49 ± 6.83 months; 34M/22F) and recruited at Beijing Children’s Hospital, Beijing, China. These infants were referred to MR imaging due to simple febrile convulsion (N = 19), convulsion (N = 21), diarrhea (N = 13), or sexual precocity (N = 3). All infants had normal neurological examinations documented in their medical records. Exclusion criteria include known nervous system diseases, or a history of neurodevelopmental, psychiatric, or systemic illness. The clinical history of each infant was carefully inspected to rule out developmental abnormalities. No brain abnormalities were detected on the MRI scans of these infants, based on the reading of an experienced pediatric radiologist (YP). All infants from this cohort-2 had PC-MRI acquired, and 48 of them had pCASL acquired. All procedures for cohort-2 were approved by Beijing Children’s Hospital Research Ethics Committee (Approval number 2016-36) and every participant’s parents provided signed consent.

Multimodal MR image acquisition

In general, multimodal brain MRIs for infants were obtained at higher resolutions in cohort-1 than in cohort-2. Brain scans of infants in cohort-1 were collected during natural sleep on a 3T Siemens Prisma scanner with a 32-channel head coil at the Children’s Hospital of Philadelphia. Brain scans of infants in cohort-2 were collected on a 3T Philips Achieva scanner with an 8-channel head coil at the Beijing Children’s Hospital. Infants from cohort-2 were under sedation with orally administered chloral hydrate at a dose of 0.5 ml/kg and no more than 10 ml in total during scanning. Previous research^53,54 suggested no significant effect of chloral hydrae on cerebral blood flow or sensory function. Earplug and headphones were used to reduce noise exposure. Below are detailed parameters for each MRI modality.

Structural MRI

Prior to perfusion MRI acquisition, structural MRIs were acquired. In cohort-1, high-resolution structural MRIs were acquired with magnetization-prepared, rapid acquisition gradient-echo T1-weighted (MPRAGE, T1w) images (TR = 2400 ms; TE = 2.24 ms; TI = 1060 ms; flip angle = 8°; sagittal FOV = 256 × 256 mm², matrix = 320 × 320, slice number = 208; resolution = isotropic 0.8 mm, scan duration = 6 min 38 s) and a variable-flip-angle turbo-spin-echo T2-weighted (Siemens SPACE, T2w) images (TR = 3200 ms; TE = 564 ms; turbo factor = 314; sagittal FOV = 256 × 240 mm², matrix = 320 × 300, slice number = 208; resolution = isotropic 0.8 mm, scan duration = 5 min 57 s). In cohort-2, structural MRI was acquired with T1-weighted images (TR = 8.28 ms; TE = 3.82 ms; TI = 1100 ms; flip angle = 12°; sagittal FOV = 200 × 200 mm², matrix = 200 × 200, slice number = 150; resolution = isotropic 1 mm, scan duration = 3.7 min).

Perfusion MRI

We collected pCASL and PC-MRI in both cohorts to measure regional and global CBF, respectively. To position the PC-MRI imaging planes, we also acquired angiogram. MR images from a representative infant were used to demonstrate pCASL, PC-MR and angiogram acquisitions (Supplementary Fig. 1a).

In cohort-1, all high-resolution pCASL scans were acquired with a 3D multi-shot, stack-of-spirals pCASL sequence^33,34 with the following parameters: TR = 4000 ms; TE = 9.4 ms; four-shot acquisition; FOV = 192 × 192 mm², matrix = 76 × 76; 48 slices; slice thickness/gap = 2.5/0 mm, effective voxel resolution = 2.53 × 2.53 × 2.5 mm³; labeling duration = 1600 ms; post labeling delay (PLD) = 1800 ms; number of controls/labels = 10 pairs; center of labeling slab located between cervical vertebrae C2 and C3 (Supplementary Fig. 1a); background suppression to suppress approximately 90% of the signal contributed from static tissue to reduce physiological noise, and scan duration = 7 min 6 s. For each pCASL scan, a pair of M₀ images were acquired and averaged for rCBF quantification with longer TR, no magnetization preparation, and the same readout scheme³³. A time-of-flight (TOF) angiogram was acquired with axial slices encompassing a slab covering the foramen magnum. The imaging parameters of the angiogram were: TR = 21 ms; TE = 3.42 ms; flip angle = 18^o, FOV = 200 × 180 mm², matrix = 256 × 144; 30 slices; slice thickness/gap = 2/0mm, effective voxel resolution = 0.78 × 0.78 × 2.0 mm³; thickness of saturation slab above the imaging slab = 40 mm, and scan duration = 47 s. The four feeding arteries, including bilateral internal carotid artery (ICA) and vertebral artery (VA), could be well visualized in the three-dimensional reconstructed angiogram on the middle and right panels of Supplementary Fig. 1a. Based on the angiogram, the imaging planes for PC-MRI of ICAs were placed at the level of the foramen magnum and the imaging planes for PC-MRI of VAs were placed between the two turns in V3 segments (at approximately the level of the C1 vertebral column) (Supplementary Fig. 1a). Image parameters of PC-MRI were: TR = 23.2 ms; TE = 5.76 ms; flip angle = 18^o; FOV = 200 × 150 mm², matrix = 256 × 120; 14 slices; slice thickness/gap = 2/0mm, effective voxel resolution = 0.78 × 0.78 × 2.0 mm³; maximum velocity encoding = 40 cm/s; non-gated; 2 repetitions, and scan duration of each artery = 33 s. As shown in the right panels of Supplementary Fig. 1a, the cross-section of the target artery with higher intensity can be observed in each PC image.

In cohort-2, infant pCASL images were acquired with a 2D multi-slice echo-planar imaging sequence with following parameters: TR = 4100 ms; TE = 15 ms; FOV = 230 × 230 mm², matrix = 84 × 84; 20 slices; slice thickness/gap = 5/0mm, effective voxel resolution = 2.74 × 2.74 × 5 mm³; labeling duration = 1650 ms; PLD = 1600 ms; number of controls/labels = 30 pairs; center of labeling slab located at the junction of spinal cord and medulla (65 mm below central slab of imaging plane), and scan duration = 4.2 min. An auxiliary scan with identical readout module to pCASL but without labeling was acquired for estimating the value of equilibrium magnetization of brain tissue. Similarly, a TOF angiogram was acquired with following parameters: TR = 20 ms; TE = 3.45 ms; flip angle = 30^o, FOV = 100 × 100 mm², matrix = 100 × 100; 30 slices; slice thickness/gap = 1/0mm, effective voxel resolution = 1 × 1 × 1 mm³, and scan duration = 28 s. PC-MRI were acquired with following parameters: TR = 20 ms; TE = 10.6 ms; flip angle = 15°; FOV = 120 × 120 mm²; matrix size = 200 × 200; single slice, effective voxel resolution = 0.6 × 0.6 × 3 mm³; maximum velocity encoding = 40 cm/s; non-gated; 4 repetitions, and scan duration of each artery = 24 s. More details for cohort-2 can be found in our previous work³⁰.

Infant MRI scan preparation for cohort-1

Infant MRI has been known to be challenging. We developed specific preparations to successfully acquire multimodal MR imaging from infants during natural sleep. Before scheduling the visit, a telephone screening is completed with the child’s caregiver to collect information, including infant’s typical nap and feeding schedule, and sleeping preferences (i.e., use of pacifier, white noise, swaddling, and crib/napping environment). The MRI scan is usually scheduled at the infant’s nap or bedtime. After scheduling, we send a confirmation email that outlines MRI visit details, what to expect, what to bring, and audio files of MRI noises to practice playing for their infant during nap time at home. Both infant and caregiver are required to change into MRI-safe hospital gowns during the imaging visit, and warmed blankets are provided for infant to be swaddled in. The MRI waiting room is used for caregivers to begin preparing for their infants to go down for their nap. Caregivers are encouraged to follow their infants’ typical nap and bedtime routines as closely as they can (i.e., diaper change, swaddle, reading a book, singing). When they are ready to enter the MRI room, the research assistant and technologist escort them into the scanning room after being screened for metal in or on the body. All lights are turned off and white noise or other soothing sounds are played inside the scanning room at the parent’s request. An MRI-compatible rocking chair is set up for the caregiver to feed and rock the infant to sleep in their arms. If rocking the infant in their arms is not possible, we allow the caregiver to let the infant fall asleep on an MRI safe bed before transitioning them to the MRI scanner table. At this time, MiniMuffs for ear protection to reduce scanner noise are placed on the infant’s ears. The caregivers are encouraged to wait 5–10 min after their infant initially falls asleep to ensure they are in a deeper sleep for transferring onto the MRI scanner table. Once asleep, the caregiver assists in gently placing the infant on the scanner table. To minimize motion during the scan, the MRI technologist positions the infant’s head in the coil using foam pads and cushions around the head.

Infant behavior and developmental functioning

Only infants from cohort-1 had their behavior and developmental functioning assessed. To examine the associations of regional CBF dynamics with measures of infant behavior, we analyzed the Bayley scales of infant and toddler development⁵⁵ that was completed during infants’ imaging visit at the Children’s Hospital of Philadelphia. This neurodevelopmental assessment was conducted by a certified neurodevelopmental psychologist (ESK), who was blinded to the infant MR findings. Fifty-three of the seventy-six infants who had high-resolution pCASL acquired in cohort-1 completed the Bayley scale to assess their behavior and developmental functioning in three domains: motor, language, and cognitive. Specifically, the cognitive scale estimates general cognitive functioning on the basis of nonverbal activities; the language scale estimates receptive communication as well as expressive communication including the ability to communicate through words and gestures; and the motor scale estimates both fine motor and gross motor. The Bayley scale is age standardized and widely used in both research and clinical settings. The three Bayley domains’ composite scores are expressed as standard scores with a mean of 100 and a standard deviation of 15, with higher scores indicating better performance. Unlike cognitive, language and motor scales reliably obtained using items administered to the infant by a certified neurodevelopmental psychologist, other two scales from Bayley (social-emotional and adaptive scales) obtained from primary caregiver heterogeneous responses to questionnaires were not included in this study.

Multimodal MR imaging processing with estimated motion

We estimated global and regional CBF from PC and pCASL perfusion MRI for infant datasets from both cohorts using the preprocessing methodologies described below.

Global CBF estimation with PC perfusion MRI

PC-MRI provides a quantitative measurement of the flow velocity of a given blood vessel. By integrating over the cross-section of the vessel, blood velocity can be converted to flow rate. The flow rates from four feeding arterials can then be used to calculate global CBF³⁵ and formulated as

$${global\; CBF}=\int {vdA}/(\rho*{brain\; volume})$$

(1)

where v is the blood flow velocity in the ICAs and Vas; A is the cross-sectional area of blood vessel with a unit of mm²; and ρ is the brain tissue density assumed as 1.06 g/mL^56,57. Brain volume was measured from structural MRI as parenchyma volume (gray matter + white matter volume).

Regional CBF estimation with pCASL perfusion MRI

After head motion correction of the pCASL perfusion MRI, we estimated rCBF using the protocol similar to that in our previous publication²⁸. Briefly, rCBF was measured using a model described in the ASL consensus paper²⁷:

$${rCBF}=\frac{6000\cdot \lambda \cdot \Delta M\cdot {e}^{\frac{{PLD}}{{T}_{1a}}}}{2\cdot \alpha \cdot {M}_{b}^{0}\cdot {T}_{1a}\cdot \left(1-{e}^{\frac{-{LabelDur}}{{T}_{1a}}}\right)}[{ml}/100g/\min ]$$

(2)

where ∆M is the dynamic-averaged signal intensity difference between in the control and label images; λ is the blood-brain partition coefficient in ml/g (0.9 ml/g⁵⁷); PLD is the post labeling delay time; LabelDur is the labeling duration; α, the labeling efficiency, is 0.86 predicted by the fitting between labeling efficiency and blood velocity in the previous study⁵⁸; T_1a is T₁ of arterial blood (1800 ms^59,60). The value of equilibrium magnetization of brain tissue (${M}_{b}^{0}$) was obtained from the M₀ images from the acquisition. The labeling efficiency α and T_1a could vary considerably across participants, especially in infants. Thus, we used global CBF from PC-MRI to calibrate rCBF measures from pCASL^28,30,58. Specifically, rCBF was calibrated by applying the scalar factor making averaged rCBF equal to global CBF. Previous validation studies have shown that this procedure enhances accuracy and reliability of rCBF quantification⁵⁸. The pCASL imaging parameters for infants, such as PLD, were selected in the interval between imaging parameters of neonates and those of children, with a PLD of 2000 ms for neonates and 1500 ms for children recommended by the ASL consensus paper²⁷. Previous studies reported the mean arterial transit time (ATT) in infants less than 1500 ms⁶⁰. Following the consensus paper²⁷ suggestion that the PLD should be larger than the ATT, a tailored PLD of 1600–1800 ms which is larger than ATT of infants was selected for the infant cohorts to minimize the effect of PLD on rCBF across subjects and brain regions.

We scanned an infant (18 months old) twice to demonstrate the reproducibility of the infant pCASL protocol in cohort-1 (see Multimodal MR imaging acquisition). To quantify the test-retest reproducibility of rCBF measurement, we further calculated the intraclass correlation coefficient (ICC) with rCBF maps estimated from these two pCASL scans (Supplementary Fig. 1b, c).

Motion metric for pCASL perfusion MRI

To quantify the in-scanner motion of pCASL, a widely used measurement of mean relative displacement (MRD) across image volumes from previous ASL and functional MRI studies e.g., refs. ^61,62 was adopted. We used the same procedures from previous studies^61,62 to obtain the MRD measurement. First, all pCASL images were motion-corrected using the MCFLIRT function from the FSL package⁶³. The re-alignment parameters, six motion parameters consisting of three translations and three rotations, estimated from MCFLIRT across pCASL volumes, were condensed to a single vector representing the root mean squared (RMS) volume-to-volume displacement of all brain voxels. This one-dimensional motion timeseries were further used to calculate the RMS displacement relative to the preceding volume (i.e., relative RMS displacement). In order to provide a summary measure of in-scanner motion for each subject, the MRD measurement was estimated as the mean value of the relative RMS displacement vector. Subjects with excessive motion (MRD > 0.5 mm) were excluded⁶¹. For each subject, any image volumes with large motion (relative RMS displacement > 0.5 mm) were removed. MRD values from all 76 infants (0.22 ± 0.13 mm) were less than 0.5 mm (Supplementary Table 3). Infant’s in-scanner motion was also included as a covariate in the subsequent analyses of developmental models of rCBF and the association analyses with infant behavior. In addition to the above-mentioned procedures that mitigate potential motion effects, high-resolution rCBF maps in cohort-1 were acquired with 3D multi-shot, stack-of-spiral pCASL with background suppression, a recommended ASL protocol by the ASL consensus paper²⁷. The protocol improves robustness against motion artifacts through stack-of-spiral acquisition for high oversampling in the center of k-space at every shot, and provides high SNR for detecting the ASL signal for rCBF calculation with the background-suppressed 3D readout^27,64.

Group-averaged rCBF maps generation

To generate the group-averaged rCBF maps across infancy, we registered rCBF maps from all infants in cohort-1 into the JHU brain template in the Montreal Neurologic Institute (MNI) space⁶⁵ (2 mm isotropic resolution with 90×109×90 matrix) utilizing the contrasts of T1w images to drive the registration. Briefly, infants’ T1w images were first co-registered to the M₀ image from pCASL perfusion MRI in the native rCBF space. Then, a 12-parameter affine registration was used to transform the co-registered T1w image of each infant to the T1w image in the template space, followed by a non-linear transformation (LDDMM: large deformation diffeomorphic metric mapping). By applying the same individual transformations, infant’s rCBF maps were then projected to the template space. The registration procedures were conducted using DiffeoMap software (www.mristudio.org). We divided the 76 infants from cohort-1 into six age groups (i.e., 0–3, 3–6, 6–9, 9–12, 12–18, and 18–28 months; Supplementary Table 3) to characterize the dynamic rCBF changes during infancy. The corresponding rCBF maps from each age group were then averaged in template space to generate the group-averaged maps. For each infant, rCBF values from the 144,237 gray matter voxels as defined by the JHU brain atlas⁶⁵ were used for subsequent analyses of voxel-wise developmental models of rCBF.

Delineate the developmental models of global and regional CBF

Identify the developmental model for age-related change of global CBF

To quantify the developmental curve of global CBF during infancy, we first identified the best model for the relationship between global CBF and age. Model selection for developmental curves of global CBF utilized the Akaike information criterion (AIC) and Bayesian information criterion (BIC) to quantify the relative goodness of fit among different candidate models^66,67, including linear, exponential, logarithmic, Poisson, and quadratic polynomial for fitting the age-related curves of global CBF during infancy (Supplementary Fig. 2). The ${AIC}=\,-2{{\mathrm{ln}}}\left(L\right)+2k$, where L = maximum likelihood of model and k = number of parameters. Since the absolute AIC value is highly specific to the dataset, the relative AIC ($\Delta {{AIC}}_{i}={{AIC}}_{i}-{{AIC}}_{\min }$) are more widely used to rank models. The value of the ΔAIC_i of the relatively best model is set to 0. The likelihood of a model can be estimated by ${e}^{-1/(2\varDelta {{AIC}}_{i})}$. The likelihood of each model divided by the sum of the likelihood for all models to estimate the AIC weight (wAIC), which can be interpreted as the probability of the model to be the best model among the candidate models (Supplementary Table 2). All analyses were conducted in R version 4.1.1⁶⁸. The AIC results showed that the logarithmic model had more than 99% probability to be the best model for characterizing global CBF development during (Supplementary Table 2). We also investigated the effects of sex and cohort on developmental curve of global CBF with a multivariate regression model (e.g., Eq. (3)).

$${global\; CBF} \sim \log \left({Age}\right)+{\beta }_{{sex}}+{\beta }_{{cohort}}$$

(3)

Identify the developmental models for age-related change of regional CBF

After analyzing the effect of age on global CBF, we sought to evaluate the spatial distribution of age effects on rCBF. Considering the logarithmic model best fitted the global CBF growth during the same developmental stage, we adopted the same model here for voxel-wise developmental models of rCBF. Infant sex and in-scanner head motion were included as covariates with each logarithmic model. Head motion was quantified as MRD (see Motion metric for pCASL perfusion MRI). This model can be summarized using the formula in Eq. (4):

$${rCBF} \sim \log \left({Age}\right)+{\beta }_{{sex}}+{\beta }_{{head\; motion}}$$

(4)

Z map was generated to demonstrate the age effects on rCBF. We used Bonferroni correction to further control the type 1 error. Given a sample of 76 subjects in cohort-1 and age effects across the 144,237 gray matter voxels, corrected P < 0.05 corresponds to uncorrected P < 3.47 × 10⁻⁷ or Z > 5.1. To visualize the dynamic distribution patterns of rCBF throughout infancy, we constructed fitted rCBF maps at specific ages ranging from 1 to 28 months with an increment of 0.5 months based on the best-fit development curves across cortex (Supplementary Figs. 6–7). In estimating the rCBF at each age, we assumed the head motion as 0.22 mm (i.e., the averaged MRD value from all subjects) and averaged the sex effect from both genders. The resulting rCBF maps were then animated to create a quantitative time-lapse movie of spatiotemporal changes of rCBF during infancy (Supplementary Movie 1).

For validation, the voxel-wise age effects on rCBF (Z map) during infancy were also replicated using the infant perfusion MRI from cohort-2. In order to evaluate the alignment of age effects on rCBF from two cohorts, Pearson correlation between two Z maps was calculated (r_oberserved). The significance of the alignment was further quantified using a nonparametric permutation test. Specifically, surrogate Z maps (N = 10,000) were generated by re-assigning voxel-wise values from Z map of cohort-2. A distribution of correlation coefficients was constructed by correlating these surrogate maps with Z map from cohort-1. The proportion of 10,000 null correlations (r_null) that exceed the observed correlation coefficient between the two Z maps (r_oberserved) is defined as an empirical P_perm value. P_perm < 0.05 indicates significant alignment of age effects on rCBF from two infant cohorts (Supplementary Fig. 5).

Segmented regression analyses of infant CBF

To identify distinct phases of CBF changes during infancy, we employed the segmented regression approach, which has also been used in previous studies e.g., refs. ^11,36 to detect break points between different phases in age-related curves of imaging measures. We used segmented regression analysis to determine the break points in the relationship between infant CBF (i.e., global and regional) and age with the Segmented (version 1.4-1) package⁶⁹ in R version 4.1.1. The approach requires the user to specify the number of breakpoints, and then identifies each break-point ___location or value. We examined a range of models with 0 to 3 break points in the segmented function for infant CBF. Increase in the number of breakpoints from 0 to 1 improved the model adjusted R² and standard error. Statistical significance was assessed using analysis of variance (ANOVA) to compare the models (e.g., linear and biphasic linear models). Increasing the number of breakpoints to 2 or 3 did not significantly improve the model. We therefore selected the biphasic model as the best fit for infant CBF in this dataset. A similar one break-point segmented regression approach was conducted for global CBF, and regional CBF from clusters identified by below clustering analyses. Segmented regression results are shown in Supplementary Table 4.

Clustering analyses of infant rCBF

To reveal the spatiotemporal heterogeneity of rCBF change in the infant brain, we parcellated brain regions using a data-driven clustering approach named non-negative matrix factorization (NMF), a reliable parcellation technique that has been widely used in neuroimaging studies e.g., refs. ^39,70. Instead of employing established cortical parcellations based on adult brains, this approach allows us to discover the infant-specific cortical topography of rCBF developmental regionalization. The NMF method generates a parts-based representation of all infants’ rCBF maps by grouping cortical voxels that change in a similar fashion, thus facilitating interpretation of rCBF developmental regionalization during infancy. Specifically, NMF factorizes a tall non-negative data matrix ${{{\bf{X}}}} \, {{\in }} \, {{{{\bf{R}}}}}^{M\times N}$ into two non-negative matrices: a component matrix ${{{\bf{W}}}} \, {{\in }} \, {{{{\bf{R}}}}}^{M\times K}$ and a coefficient matrix ${{{\bf{H}}}} \, {{\in }} \, {{{{\bf{R}}}}}^{K\times N}$. In other words, the new variables ${{{\bf{W}}}}$, termed hereafter components or clusters, were constructed as a linear, non-negative combination of the original variables. NMF aggregated variance in the components by positively weighting original variables (i.e., rCBF) that tend to co-vary. The scalar K represents the number of components (i.e., clusters), which is usually small (i.e., $K\ll M$ and $K\ll N$). In contrast to clustering approaches that assign each cortical voxel to a single component, NMF yields a soft (probabilistic) parcellation. This soft parcellation can be converted into hard (discrete) regionalization by assigning each cortical voxel to a specific component according to its highest weight. Estimating the components with NMF can be formulated as:

$${\min }_{{{{\bf{W}}}}}{|\left|{{{\bf{X}}}}-{{{\bf{WH}}}}\right||}_{F}^{2} \, {s}.t.\,{{{\bf{W}}}},\,{{{\bf{H}}}}\ge 0,\,{{{{\bf{W}}}}}^{{{{\rm{T}}}}}{{{\bf{W}}}}={{{\bf{I}}}},\,{{{\bf{H}}}}={{{{\bf{W}}}}}^{{{{\rm{T}}}}}{{{\bf{X}}}}$$

(5)

In this study, ${{{\bf{X}}}}{{\in }}{{{{\bf{R}}}}}^{M\times N}$ is a large non-negative data matrix containing rCBF values from all infants registered to the template space, where M and N represent the number of cortical voxels and subjects, respectively. To find an appropriate number of components (i.e., K) for rCBF growth, we adopted Davies-Bouldin criterion⁷¹ and Calinski-Harabasz criterion⁷² to evaluate the performance of clustering analyses. The Davies-Bouldin index represents average similarity measure of each cluster to its most similar cluster, and a lower value of the index means the clusters are well separated. The Calinski-Harabasz index represents the ratio of the sum between-cluster dispersion and within-cluster dispersion for all clusters, and higher value of the index means the clusters are dense and well separated. We selected component number of 3 with the best criteria scores for categorizing rCBF development, after testing the clustering analyses with K ranging from 2 to 20 (Supplementary Fig. 8). The components can be represented as soft cluster assignment maps (Fig. 4a with K = 3), where the distribution of the coefficient values of rCBF exhibiting similar statistical properties peaks for the same component.

Associations of infant rCBF dynamics with infant behavior

We evaluated associations between infant rCBF and developmental functioning at each cortical voxel using generalized additive models (GAMs) with penalized splines to account for linear and nonlinear effect of age⁴⁰. Infant sex, in-scanner head motion, and family socioeconomic status (SES) were included as covariates in the model (Eq. (6)). Here, infant developmental functioning was quantified with Bayley scale across three domains: motor, language, and cognitive (see Infant behavior and developmental functioning for details). Family SES was derived from parental education and occupation using the Amherst modification of the Hollingshead two-factor index⁷³. Among the infants who had high-resolution pCASL acquired in cohort-1, only 49 of them had both Bayley scale and family SES information completed and data from these 49 infants were included in the association analyses. Cohort-2 was not included in the rCBF-behavior analyses. We used the mgcv (version 1.8-39) package in R version 4.1.1 for the analyses.

$${rCBF} \sim {Behavior}+s\left({Age}\right)+{\beta }_{{sex}}+{\beta }_{{head\; motion}}+{\beta }_{{SES}}$$

(6)

To avoid false positive results, the threshold for voxel-wise across-subject correlation maps was set to a corrected P < 0.05, which corresponded to an uncorrected single voxel significance level of P < 0.05 (t > 2.02) and a minimum cluster extent of 100 continuous voxels (800 mm³). The cortical region to which each cluster belongs was identified with infant-specific NMF components (see Clustering analyses of infant rCBF for details) and JHU brain atlas⁶³ (Supplementary Table 5). The rCBF values of these cluster voxels were measured to calculate the average rCBF at the identified clusters for each infant. Subsequently, cluster-level across-subject association analyses were performed on the averaged rCBF values and 3 Bayley scores using Eq. (6). Familywise error-corrected significance threshold at cluster level was set at 0.05/3 = 0.0167. Association between cluster-level averaged rCBF and infant developmental functioning was also conducted while controlling for covariates (i.e., age, sex, in-scanner head motion, and family SES).

Associations of infant rCBF topography with cerebral metabolism in adult brain

To assess the spatial alignment between rCBF and cerebral metabolism, we conducted a Pearson correlation analysis between the rCBF maps and the cerebral metabolism map in the adult brain. To assess the developmental reorganization of infant rCBF, we further conducted Pearson correlation analyses between the fitted rCBF maps across infancy (see Identify the developmental models for age-related change of regional CBF for details) and the cerebral metabolism map in adult brain. Cerebral metabolism rate of glucose consumption (CMRglc) with a unit of µmol/100 g/min, estimated by ¹⁸F-flurodeoxyglucose; fluorodeoxyglucose-positron emission tomography (FDG-PET), was used to quantify the cerebral metabolism. The group-averaged (N = 28) CMRglc map for adult brain was obtained from the publicly available dataset (see details in ref. ⁷⁴; https://github.com/eshoko/COMET) and was registered to the template of infant rCBF maps using linear and nonlinear registrations from FSL package. Association analyses were then performed between the fitted rCBF maps and registered CMRglc map. Permutation tests with 10,000 surrogate maps were used to estimate the significance of these alignments across infancy. Details on permutation testing can be found in the above section on developmental models for age-related change of regional CBF, which describes the alignment of age effects on rCBF across two cohorts.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Data availability

All data required for reproducing our findings are publicly available, including the global cerebral blood flow measures and regional cerebral blood flow maps, infant behavioral scores and the data for visualizing main figures. They are available and can be downloaded via the GitHub repository (https://github.com/haohuanglab/infant_rCBF_dynamics). The complete infant behavioral tests datasets used in this study are available in the NIMH Data Archive (NDA) (collection C2508) at https://nda.nih.gov/edit_collection.html?id=2508. The adult cerebral metabolism rate of glucose consumption data used in this study is available at https://github.com/eshoko/COMET. Source data are provided with this paper.

Code availability

Software packages used in this study include DiffeoMap (https://www.mristudio.org), FSL package (v5.0.9) (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki), Matlab (R2017a; MathWorks, Inc), non-negative matrix factorization package (https://github.com/asotiras/brainparts), R 4.1.1 (https://www.r-project.org), Segmented (version 1.4-1) package (https://cran.r-project.org/web/packages/segmented/index.html), and mgcv (version 1.8-39) package (https://cran.r-project.org/web/packages/mgcv/index.html). The code used in the analysis and related documentation is available at the GitHub repository (https://github.com/haohuanglab/infant_rCBF_dynamics) and is also linked to Zenodo⁷⁵.

References

Bethlehem, R. A. I. et al. Brain charts for the human lifespan. Nature 604, 525–533 (2022).
Article PubMed PubMed Central CAS Google Scholar
Ouyang, M., Dubois, J., Yu, Q., Mukherjee, P. & Huang, H. Delineation of early brain development from fetuses to infants with diffusion MRI and beyond. Neuroimage 185, 836–850 (2019).
Article PubMed Google Scholar
Attwell, D. et al. Glial and neuronal control of brain blood flow. Nature 468, 232–243 (2010).
Article ADS PubMed PubMed Central CAS Google Scholar
Kennedy, C. & Sokoloff, L. An adaptation of the nitrous oxide method to the study of the cerebral circulation in children; normal values for cerebral blood flow and cerebral metabolic rate in childhood. J. Clin. Investig. 36, 1130–1137 (1957).
Article PubMed PubMed Central CAS Google Scholar
Kuzawa, C. W. et al. Metabolic costs and evolutionary implications of human brain development. Proc. Natl Acad. Sci. USA 111, 13010–13015 (2014).
Article ADS PubMed PubMed Central CAS Google Scholar
Silbereis, J. C., Pochareddy, S., Zhu, Y., Li, M. & Sestan, N. The cellular and molecular landscapes of the developing human central nervous system. Neuron 89, 248–268 (2016).
Article PubMed PubMed Central CAS Google Scholar
Vaishnavi, S. N. et al. Regional aerobic glycolysis in the human brain. Proc. Natl Acad. Sci. USA 107, 17757–17762 (2010).
Article ADS PubMed PubMed Central CAS Google Scholar
Gur, R. C. et al. Regional differences in the coupling between resting cerebral blood flow and metabolism may indicate action preparedness as a default state. Cereb. Cortex 19, 375–382 (2009).
Article PubMed Google Scholar
Dobbing, J. & Sands, J. Comparative aspects of the brain growth spurt. Early Hum. Dev. 3, 79–83 (1979).
Article PubMed CAS Google Scholar
Volpe, J. J. Brain injury in premature infants: a complex amalgam of destructive and developmental disturbances. Lancet Neurol. 8, 110–124 (2009).
Article PubMed PubMed Central Google Scholar
Ouyang, M. et al. Differential cortical microstructural maturation in the preterm human brain with diffusion kurtosis and tensor imaging. Proc. Natl Acad. Sci. USA 116, 4681–4688 (2019).
Article ADS PubMed PubMed Central CAS Google Scholar
Yakovlev, P. I. The myelogenetic cycles of regional maturation of the brain. Regional development of the brain in early life, 3–70 (1967).
Sydnor, V. J. et al. Neurodevelopment of the association cortices: patterns, mechanisms, and implications for psychopathology. Neuron 109, 2820–2846 (2021).
Article PubMed PubMed Central CAS Google Scholar
Huang, H. Imaging the Infant Brain. (Oxford Research Encyclopedia of Psychology, 2022).
Hill, J. et al. Similar patterns of cortical expansion during human development and evolution. Proc. Natl Acad. Sci. USA 107, 13135–13140 (2010).
Article ADS PubMed PubMed Central CAS Google Scholar
Wang, F. et al. Developmental topography of cortical thickness during infancy. Proc. Natl Acad. Sci. USA 116, 15855–15860 (2019).
Article ADS PubMed PubMed Central CAS Google Scholar
Huang, H. et al. Anatomical characterization of human fetal brain development with diffusion tensor magnetic resonance imaging. J. Neurosci. 29, 4263–4273 (2009).
Article PubMed PubMed Central CAS Google Scholar
Huang, H. et al. White and gray matter development in human fetal, newborn and pediatric brains. Neuroimage 33, 27–38 (2006).
Article PubMed Google Scholar
Vasung, L. et al. Development of axonal pathways in the human fetal fronto-limbic brain: histochemical characterization and diffusion tensor imaging. J. Anat. 217, 400–417 (2010).
Article PubMed PubMed Central Google Scholar
Cao, M. et al. Early development of functional network segregation revealed by connectomic analysis of the preterm human brain. Cereb. Cortex 27, 1949–1963 (2017).
PubMed Google Scholar
Smyser, C. D. et al. Longitudinal analysis of neural network development in preterm infants. Cereb. Cortex 20, 2852–2862 (2010).
Article PubMed PubMed Central Google Scholar
Gao, W. et al. Evidence on the emergence of the brain’s default network from 2-week-old to 2-year-old healthy pediatric subjects. Proc. Natl Acad. Sci. USA 106, 6790–6795 (2009).
Article ADS PubMed PubMed Central CAS Google Scholar
Chugani, H. T. & Phelps, M. E. Maturational changes in cerebral function in infants determined by 18FDG positron emission tomography. Science 231, 840–843 (1986).
Article ADS PubMed CAS Google Scholar
Chugani, H. T., Phelps, M. E. & Mazziotta, J. C. Positron emission tomography study of human brain functional development. Ann. Neurol. 22, 487–497 (1987).
Article PubMed CAS Google Scholar
Drew, P. J., Mateo, C., Turner, K. L., Yu, X. & Kleinfeld, D. Ultra-slow oscillations in fMRI and resting-state connectivity: neuronal and vascular contributions and technical confounds. Neuron 107, 782–804 (2020).
Article PubMed PubMed Central CAS Google Scholar
Liu, T. T. Noise contributions to the fMRI signal: an overview. Neuroimage 143, 141–151 (2016).
Article PubMed Google Scholar
Alsop, D. C. et al. Recommended implementation of arterial spin-labeled perfusion MRI for clinical applications: A consensus of the ISMRM perfusion study group and the European consortium for ASL in dementia. Magn. Reson Med. 73, 102–116 (2015).
Article ADS PubMed Google Scholar
Ouyang, M. et al. Heterogeneous increases of regional cerebral blood flow during preterm brain development: Preliminary assessment with pseudo-continuous arterial spin labeled perfusion MRI. Neuroimage 147, 233–242 (2017).
Article PubMed Google Scholar
Wang, Z. et al. Assessment of functional development in normal infant brain using arterial spin labeled perfusion MRI. Neuroimage 39, 973–978 (2008).
Article PubMed Google Scholar
Yu, Q. et al. Infant brain regional cerebral blood flow increases supporting emergence of the default-mode network. Elife 12, e78397 (2023).
Lemaitre, H. et al. Rest functional brain maturation during the first year of life. Cereb. Cortex 31, 1776–1785 (2021).
Article PubMed Google Scholar
Jain, V. et al. Longitudinal reproducibility and accuracy of pseudo-continuous arterial spin-labeled perfusion MR imaging in typically developing children. Radiology 263, 527–536 (2012).
Article PubMed PubMed Central Google Scholar
Chang, Y. V., Vidorreta, M., Wang, Z. & Detre, J. A. 3D-accelerated, stack-of-spirals acquisitions and reconstruction of arterial spin labeling MRI. Magn. Reson Med. 78, 1405–1419 (2017).
Article PubMed Google Scholar
Vidorreta, M. et al. Whole-brain background-suppressed pCASL MRI with 1D-accelerated 3D RARE Stack-Of-Spirals readout. PLoS One 12, e0183762 (2017).
Article PubMed PubMed Central Google Scholar
Spilt, A. et al. Reproducibility of total cerebral blood flow measurements using phase contrast magnetic resonance imaging. J. Magn. Reson Imaging 16, 1–5 (2002).
Article PubMed Google Scholar
Pontzer, H. et al. Daily energy expenditure through the human life course. Science 373, 808–812 (2021).
Article ADS PubMed PubMed Central CAS Google Scholar
Chen, Y., Wang, D. J. & Detre, J. A. Test-retest reliability of arterial spin labeling with common labeling strategies. J. Magn. Reson Imaging 33, 940–949 (2011).
Article PubMed PubMed Central Google Scholar
Yeo, B. T. et al. The organization of the human cerebral cortex estimated by intrinsic functional connectivity. J. Neurophysiol. 106, 1125–1165 (2011).
Article PubMed Google Scholar
Lee, D. D. & Seung, H. S. Learning the parts of objects by non-negative matrix factorization. Nature 401, 788–791 (1999).
Article ADS PubMed CAS Google Scholar
Wood, S. N. Stable and efficient multiple smoothing parameter estimation for generalized additive models. J. Am. Stat. Assoc. 99, 673–686 (2004).
Article MathSciNet Google Scholar
Moskowitz, M. A., Lo, E. H. & Iadecola, C. The science of stroke: mechanisms in search of treatments. Neuron 67, 181–198 (2010).
Article PubMed PubMed Central CAS Google Scholar
Huttenlocher, P. R. & Dabholkar, A. S. Regional differences in synaptogenesis in human cerebral cortex. J. Comp. Neurol. 387, 167–178 (1997).
Article PubMed CAS Google Scholar
Lebenberg, J. et al. Mapping the asynchrony of cortical maturation in the infant brain: A MRI multi-parametric clustering approach. Neuroimage 185, 641–653 (2019).
Article PubMed CAS Google Scholar
Zhu, Y. et al. Spatiotemporal transcriptomic divergence across human and macaque brain development. Science 362, eaat8077 (2018).
Li, M. et al. Integrative functional genomic analysis of human brain development and neuropsychiatric risks. Science 362, eaat7615 (2018).
Attwell, D. & Laughlin, S. B. An energy budget for signaling in the grey matter of the brain. J. Cereb. Blood Flow. Metab. 21, 1133–1145 (2001).
Article PubMed CAS Google Scholar
Liang, X., Zou, Q., He, Y. & Yang, Y. Coupling of functional connectivity and regional cerebral blood flow reveals a physiological basis for network hubs of the human brain. Proc. Natl Acad. Sci. USA 110, 1929–1934 (2013).
Article ADS PubMed PubMed Central CAS Google Scholar
Choi, D., Yeung, H. H. & Werker, J. F. Sensorimotor fundations of speech perception in infancy. Trends Cogn. Sci. 27, 773–784 (2023).
Article PubMed Google Scholar
van der Velde, B., White, T. & Kemner, C. The emergence of a theta social brain network during infancy. Neuroimage 240, 118298 (2021).
Article PubMed Google Scholar
Bayley, N., & Aylward, G. P. Bayley Scales of Infant and Toddler Development: Technical manual (4th ed.). (Bloomington, MN: NCS Pearson, 2019).
Estes, A. et al. Behavioral, cognitive, and adaptive development in infants with autism spectrum disorder in the first 2 years of life. J. Neurodev. Disord. 7, 24 (2015).
Article PubMed PubMed Central Google Scholar
Lindner, T. et al. Current state and guidance on arterial spin labeling perfusion MRI in clinical neuroimaging. Magn. Reson Med. 89, 2024–2047 (2023).
Article PubMed PubMed Central Google Scholar
Li, A., Gong, L. & Xu, F. Brain-state–independent neural representation of peripheral stimulation in rat olfactory bulb. Proc. Natl Acad. Sci. USA 108, 5087–5092 (2011).
Article ADS PubMed PubMed Central CAS Google Scholar
Suzuki, C., Kosugi, M. & Magata, Y. Conscious rat PET imaging with soft immobilization for quantitation of brain functions: comprehensive assessment of anesthesia effects on cerebral blood flow and metabolism. EJNMMI Res. 11, 1–11 (2021).
Article Google Scholar
Bayley, N. Bayley scales of infant development. (TX: Harcourt Assessment, San Antonio, 2006).
Dittmer, J. C. & Dawson, R. M. The isolation of a new lipid, triphosphoinositide, and monophosphoinositide from ox brain. Biochem. J. 81, 535–540 (1961).
Article PubMed PubMed Central CAS Google Scholar
Herscovitch, P. & Raichle, M. E. What is the correct value for the brain–blood partition coefficient for water? J. Cereb. Blood Flow. Metab. 5, 65–69 (1985).
Article PubMed CAS Google Scholar
Aslan, S. et al. Estimation of labeling efficiency in pseudocontinuous arterial spin labeling. Magn. Reson Med. 63, 765–771 (2010).
Article PubMed PubMed Central Google Scholar
Liu, P. et al. T1 and T2 values of human neonatal blood at 3 Tesla: dependence on hematocrit, oxygenation, and temperature. Magn. Reson Med. 75, 1730–1735 (2016).
Article PubMed Google Scholar
Varela, M. et al. A method for rapid in vivo measurement of blood T1. NMR Biomed. 24, 80–88 (2011).
Article PubMed Google Scholar
Satterthwaite, T. D. et al. Impact of puberty on the evolution of cerebral perfusion during adolescence. Proc. Natl Acad. Sci. USA 111, 8643–8648 (2014).
Article ADS PubMed PubMed Central CAS Google Scholar
Satterthwaite, T. D. et al. Heterogeneous impact of motion on fundamental patterns of developmental changes in functional connectivity during youth. Neuroimage 83, 45–57 (2013).
Article PubMed Google Scholar
Jenkinson, M., Bannister, P., Brady, M. & Smith, S. Improved optimization for the robust and accurate linear registration and motion correction of brain images. Neuroimage 17, 825–841 (2002).
Article PubMed Google Scholar
Vidorreta, M. et al. Comparison of 2D and 3D single-short ASL perfusion fMRI seuqences. Neuroimage 66, 662–671 (2012).
Article PubMed Google Scholar
Oishi, K. et al. Atlas-based whole brain white matter analysis using large deformation diffeomorphic metric mapping: application to normal elderly and Alzheimer’s disease participants. Neuroimage 46, 486–499 (2009).
Article PubMed Google Scholar
Akaike, H. Information theory and an extension of the maximum likelihood principle. Selected papers of hirotugu akaike, 199–213 (Springer, 1998).
Burnham, K. P. & Anderson, D. R. Multimodel inference: understanding AIC and BIC in model selection. Sociological methods Res. 33, 261–304 (2004).
Article MathSciNet Google Scholar
Team, R. C. R: A language and environment for statistical computing. R Foundation for Statistical Computing. (Vienna, Austria, 2020).
Muggeo, V. M. Segmented: an R package to fit regression models with broken-line relationships. R. N. 8, 20–25 (2008).
Google Scholar
Sotiras, A., Resnick, S. M. & Davatzikos, C. Finding imaging patterns of structural covariance via non-negative matrix factorization. Neuroimage 108, 1–16 (2015).
Article PubMed Google Scholar
Davies, D. L. & Bouldin, D. W. A cluster separation measure. IEEE Trans. Pattern Anal. Mach. Intell. PAMI-1, 224–227 (1979).
Article Google Scholar
Caliński, T. & Harabasz, J. A dendrite method for cluster analysis. Commun. Stat. Theory Methods 3, 1–27 (1974).
Article MathSciNet Google Scholar
Hollingshead A. B. Two factor index of social position. (New Haven: Yale University Press, 1957).
Shokri-Kojori, E. et al. Correspondence between cerebral glucose metabolism and BOLD reveals relative power and cost in human brain. Nat. Commun. 10, 690 (2019).
Article ADS PubMed PubMed Central CAS Google Scholar
Ouyang, M. et al. Spatiotemporal cerebral blood flow dynamics underlies emergence of the limbic-sensorimotor-association cortical gradient in human infancy. Zenodo, https://doi.org/10.5281/zenodo.13625828 (2024).

Download references

Acknowledgements

This work was supported by grants from National Institute of Health: R01MH092535 (H.H.), R01MH125333 (H.H.), R01EB031284 (H.H.), R01MH129981 (H.H.), R21MH123930 (M.O.), R01EB031080 (J.A.D.), R01HD093776 (J.C.E.), R01MH107506 (J.C.E.), and P50HD105354 (H.H., J.C.E.).

Author information

Authors and Affiliations

Department of Radiology, Children’s Hospital of Philadelphia, Philadelphia, PA, USA
Minhui Ouyang, Jessica L. Hyland, Kay L. Sindabizera, Emily S. Kuschner, J. Christopher Edgar & Hao Huang
Department of Radiology, Perelman School of Medicine, University of Pennsylvania, Philadelphia, PA, USA
Minhui Ouyang, John A. Detre, J. Christopher Edgar & Hao Huang
Department of Neurology, Perelman School of Medicine, University of Pennsylvania, Philadelphia, PA, USA
John A. Detre
Department of Psychiatry, Perelman School of Medicine, University of Pennsylvania, Philadelphia, PA, USA
Emily S. Kuschner
Department of Radiology, Beijing Children’s Hospital, Capital Medical University, Beijing, China
Yun Peng

Authors

Minhui Ouyang
View author publications
Search author on:PubMed Google Scholar
John A. Detre
View author publications
Search author on:PubMed Google Scholar
Jessica L. Hyland
View author publications
Search author on:PubMed Google Scholar
Kay L. Sindabizera
View author publications
Search author on:PubMed Google Scholar
Emily S. Kuschner
View author publications
Search author on:PubMed Google Scholar
J. Christopher Edgar
View author publications
Search author on:PubMed Google Scholar
Yun Peng
View author publications
Search author on:PubMed Google Scholar
Hao Huang
View author publications
Search author on:PubMed Google Scholar

Contributions

M.O. and H.H. designed the research. M.O., J.A.D., J.L.H, K.L.S., E.S.K., J.C.E., Y.P., and H.H. performed the research. M.O., J.A.D. and H.H. contributed new methodology and analytic tools. M.O. and H.H. supervised the collection of neuroimaging data. M.O. analyzed data, and generated figures and tables. M.O. wrote the original draft and all authors (M.O., J.A.D., J.L.H, K.L.S., E.S.K., J.C.E., Y.P., and H.H.) reviewed and revised the final darft.

Corresponding author

Correspondence to Hao Huang.

Ethics declarations

Competing interests

The authors declare no competing interests.

Peer review

Peer review information

Nature Communications thanks the anonymous reviewer(s) for their contribution to the peer review of this work. A peer review file is available.

Additional information

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Supplementary Information

Peer Review File

Description of Additional Supplementary Files

Supplementary Move 1: 4D spatiotemporal dynamics of infant regional cerebral blood flow from 1 to 28 months

Reporting Summary

Source data

Source Data

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Cite this article

Ouyang, M., Detre, J.A., Hyland, J.L. et al. Spatiotemporal cerebral blood flow dynamics underlies emergence of the limbic-sensorimotor-association cortical gradient in human infancy. Nat Commun 15, 8944 (2024). https://doi.org/10.1038/s41467-024-53354-7

Download citation

Received: 24 March 2024
Accepted: 02 October 2024
Published: 17 October 2024
DOI: https://doi.org/10.1038/s41467-024-53354-7

This article is cited by

Advances in magnetic resonance imaging of the developing brain and its applications in pediatrics
- Rui-Ke Chen
- Ming-Yang Li
- Dan Wu
World Journal of Pediatrics (2025)