Characterizing genetic intra-tumor heterogeneity across 2,658 human cancer genomes.

Intra-tumor heterogeneity (ITH) is a mechanism of therapeutic resistance and therefore an important clinical challenge. However, the extent, origin, and drivers of ITH across cancer types are poorly understood. To address this, we extensively characterize ITH across whole-genome sequences of 2,658 cancer samples spanning 38 cancer types. Nearly all informative samples (95.1%) contain evidence of distinct subclonal expansions with frequent branching relationships between subclones. We observe positive selection of subclonal driver mutations across most cancer types and identify cancer type-specific subclonal patterns of driver gene mutations, fusions, structural variants, and copy number alterations as well as dynamic changes in mutational processes between subclonal expansions. Our results underline the importance of ITH and its drivers in tumor evolution and provide a pan-cancer resource of comprehensively annotated subclonal events from whole-genome sequencing data.

Identification and validation of supervariants reveal novel loci associated with human white matter microstructure.

As an essential part of the central nervous system, white matter coordinates communications between different brain regions and is related to a wide range of neurodegenerative and neuropsychiatric disorders. Previous genome-wide association studies (GWASs) have uncovered loci associated with white matter microstructure. However, GWASs suffer from limited reproducibility and difficulties in detecting multi-single-nucleotide polymorphism (multi-SNP) and epistatic effects. In this study, we adopt the concept of supervariants, a combination of alleles in multiple loci, to account for potential multi-SNP effects. We perform supervariant identification and validation to identify loci associated with 22 white matter fractional anisotropy phenotypes derived from diffusion tensor imaging. To increase reproducibility, we use United Kingdom (UK) Biobank White British (n = 30,842) data for discovery and internal validation, and UK Biobank White but non-British (n = 1927) data, Europeans from the Adolescent Brain Cognitive Development study (n = 4399) data, and Europeans from the Human Connectome Project (n = 319) data for external validation. We identify 23 novel loci on the discovery set that have not been reported in the previous GWASs on white matter microstructure. Among them, three supervariants on genomic regions 5q35.1, 8p21.2, and 19q13.32 have P-values lower than 0.05 in the meta-analysis of the three independent validation data sets. These supervariants contain genetic variants located in genes that have been related to brain structures, cognitive functions, and neuropsychiatric diseases. Our findings provide a better understanding of the genetic architecture underlying white matter microstructure.

A comprehensive patient-specific prediction model for temporomandibular joint osteoarthritis progression.

Temporomandibular joint osteoarthritis (TMJ OA) is a prevalent degenerative disease characterized by chronic pain and impaired jaw function. The complexity of TMJ OA has hindered the development of prognostic tools, posing a significant challenge in timely, patient-specific management. Addressing this gap, our research employs a comprehensive, multidimensional approach to advance TMJ OA prognostication. We conducted a prospective study with 106 subjects, 74 of whom were followed up after 2 to 3 y of conservative treatment. Central to our methodology is the development of an innovative, open-source predictive modeling framework, the Ensemble via Hierarchical Predictions through Nested cross-validation tool (EHPN). This framework synergistically integrates 18 feature selection, statistical, and machine learning methods to yield an accuracy of 0.87, with an area under the ROC curve of 0.72 and an F1 score of 0.82. Our study, beyond technical advancements, emphasizes the global impact of TMJ OA, recognizing its unique demographic occurrence. We highlight key factors influencing TMJ OA progression. Using SHAP analysis, we identified personalized prognostic predictors: lower values of headache, lower back pain, restless sleep, condyle high gray level-GL-run emphasis, articular fossa GL nonuniformity, and long-run low GL emphasis; and higher values of superior joint space, mouth opening, saliva Vascular-endothelium-growth-factor, Matrix-metalloproteinase-7, serum Epithelial-neutrophil-activating-peptide, and age indicate recovery likelihood. Our multidimensional and multimodal EHPN tool enhances clinicians' decision-making, offering a transformative translational infrastructure. The EHPN model stands as a significant contribution to precision medicine, offering a paradigm shift in the management of temporomandibular disorders and potentially influencing broader applications in personalized healthcare.

The evolutionary history of 2,658 cancers.

Cancer develops through a process of somatic evolution1,2. Sequencing data from a single biopsy represent a snapshot of this process that can reveal the timing of specific genomic aberrations and the changing influence of mutational processes3. Here, by whole-genome sequencing analysis of 2,658 cancers as part of the Pan-Cancer Analysis of Whole Genomes (PCAWG) Consortium of the International Cancer Genome Consortium (ICGC) and The Cancer Genome Atlas (TCGA)4, we reconstruct the life history and evolution of mutational processes and driver mutation sequences of 38 types of cancer. Early oncogenesis is characterized by mutations in a constrained set of driver genes, and specific copy number gains, such as trisomy 7 in glioblastoma and isochromosome 17q in medulloblastoma. The mutational spectrum changes significantly throughout tumour evolution in 40% of samples. A nearly fourfold diversification of driver genes and increased genomic instability are features of later stages. Copy number alterations often occur in mitotic crises, and lead to simultaneous gains of chromosomal segments. Timing analyses suggest that driver mutations often precede diagnosis by many years, if not decades. Together, these results determine the evolutionary trajectories of cancer, and highlight opportunities for early cancer detection.

FPLS-DC: functional partial least squares through distance covariance for imaging genetics.

MOTIVATION: Imaging genetics integrates imaging and genetic techniques to examine how genetic variations influence the function and structure of organs like the brain or heart, providing insights into their impact on behavior and disease phenotypes. The use of organ-wide imaging endophenotypes has increasingly been used to identify potential genes associated with complex disorders. However, analyzing organ-wide imaging data alongside genetic data presents two significant challenges: high dimensionality and complex relationships. To address these challenges, we propose a novel, nonlinear inference framework designed to partially mitigate these issues. RESULTS: We propose a functional partial least squares through distance covariance (FPLS-DC) framework for efficient genome wide analyses of imaging phenotypes. It consists of two components. The first component utilizes the FPLS-derived base functions to reduce image dimensionality while screening genetic markers. The second component maximizes the distance correlation between genetic markers and projected imaging data, which is a linear combination of the FPLS-basis functions, using simulated annealing algorithm. In addition, we proposed an iterative FPLS-DC method based on FPLS-DC framework, which effectively overcomes the influence of inter-gene correlation on inference analysis. We efficiently approximate the null distribution of test statistics using a gamma approximation. Compared to existing methods, FPLS-DC offers computational and statistical efficiency for handling large-scale imaging genetics. In real-world applications, our method successfully detected genetic variants associated with the hippocampus, demonstrating its value as a statistical toolbox for imaging genetic studies. AVAILABILITY AND IMPLEMENTATION: The FPLS-DC method we propose opens up new research avenues and offers valuable insights for analyzing functional and high-dimensional data. In addition, it serves as a useful tool for scientific analysis in practical applications within the field of imaging genetics research. The R package FPLS-DC is available in Github: https://github.com/BIG-S2/FPLSDC.

The interaction effects of age, APOE and common environmental risk factors on human brain structure.

Mounting evidence suggests considerable diversity in brain aging trajectories, primarily arising from the complex interplay between age, genetic, and environmental risk factors, leading to distinct patterns of micro- and macro-cerebral aging. The underlying mechanisms of such effects still remain unclear. We conducted a comprehensive association analysis between cerebral structural measures and prevalent risk factors, using data from 36,969 UK Biobank subjects aged 44-81. Participants were assessed for brain volume, white matter diffusivity, Apolipoprotein E (APOE) genotypes, polygenic risk scores, lifestyles, and socioeconomic status. We examined genetic and environmental effects and their interactions with age and sex, and identified 726 signals, with education, alcohol, and smoking affecting most brain regions. Our analysis revealed negative age-APOE-ε4 and positive age-APOE-ε2 interaction effects, respectively, especially in females on the volume of amygdala, positive age-sex-APOE-ε4 interaction on the cerebellar volume, positive age-excessive-alcohol interaction effect on the mean diffusivity of the splenium of the corpus callosum, positive age-healthy-diet interaction effect on the paracentral volume, and negative APOE-ε4-moderate-alcohol interaction effects on the axial diffusivity of the superior fronto-occipital fasciculus. These findings highlight the need of considering age, sex, genetic, and environmental joint effects in elucidating normal or abnormal brain aging.

Mapping Genetic Topography of Cortical Thickness and Surface Area in Neonatal Brains.

Adult twin neuroimaging studies have revealed that cortical thickness (CT) and surface area (SA) are differentially influenced by genetic information, leading to their spatially distinct genetic patterning and topography. However, the postnatal origins of the genetic topography of CT and SA remain unclear, given the dramatic cortical development from neonates to adults. To fill this critical gap, this study unprecedentedly explored how genetic information differentially regulates the spatial topography of CT and SA in the neonatal brain by leveraging brain magnetic resonance (MR) images from 202 twin neonates with minimal influence by the complicated postnatal environmental factors. We capitalized on infant-dedicated computational tools and a data-driven spectral clustering method to parcellate the cerebral cortex into a set of distinct regions purely according to the genetic correlation of cortical vertices in terms of CT and SA, respectively, and accordingly created the first genetically informed cortical parcellation maps of neonatal brains. Both genetic parcellation maps exhibit bilaterally symmetric and hierarchical patterns, but distinct spatial layouts. For CT, regions with closer genetic relationships demonstrate an anterior-posterior (A-P) division, while for SA, regions with greater genetic proximity are typically within the same lobe. Certain genetically informed regions exhibit strong similarities between neonates and adults, with the most striking similarities in the medial surface in terms of SA, despite their overall substantial differences in genetic parcellation maps. These results greatly advance our understanding of the development of genetic influences on the spatial patterning of cortical morphology.SIGNIFICANCE STATEMENT Genetic influences on cortical thickness (CT) and surface area (SA) are complex and could evolve throughout the lifespan. However, studies revealing distinct genetic topography of CT and SA have been limited to adults. Using brain structural magnetic resonance (MR) images of twins, we unprecedentedly discovered the distinct genetically-informed parcellation maps of CT and SA in neonatal brains, respectively. Each genetic parcellation map comprises a distinct spatial layout of cortical regions, where vertices within the same region share high genetic correlation. These genetic parcellation maps of CT and SA of neonates largely differ from those of adults, despite their highly remarkable similarities in the medial cortex of SA. These discoveries provide important insights into the genetic organization of the early cerebral cortex development.

The sequence kernel association test for multicategorical outcomes.

Disease heterogeneity is ubiquitous in biomedical and clinical studies. In genetic studies, researchers are increasingly interested in understanding the distinct genetic underpinning of subtypes of diseases. However, existing set-based analysis methods for genome-wide association studies are either inadequate or inefficient to handle such multicategorical outcomes. In this paper, we proposed a novel set-based association analysis method, sequence kernel association test (SKAT)-MC, the sequence kernel association test for multicategorical outcomes (nominal or ordinal), which jointly evaluates the relationship between a set of variants (common and rare) and disease subtypes. Through comprehensive simulation studies, we showed that SKAT-MC effectively preserves the nominal type I error rate while substantially increases the statistical power compared to existing methods under various scenarios. We applied SKAT-MC to the Polish breast cancer study (PBCS), and identified gene FGFR2 was significantly associated with estrogen receptor (ER)+ and ER- breast cancer subtypes. We also investigated educational attainment using UK Biobank data ( N = 127 , 127 $N=127,127$ ) with SKAT-MC, and identified 21 significant genes in the genome. Consequently, SKAT-MC is a powerful and efficient analysis tool for genetic association studies with multicategorical outcomes. A freely distributed R package SKAT-MC can be accessed at https://github.com/Zhiwen-Owen-Jiang/SKATMC.

Estimation for the bivariate quantile varying coefficient model with application to diffusion tensor imaging data analysis.

Despite interest in the joint modeling of multiple functional responses such as diffusion properties in neuroimaging, robust statistical methods appropriate for this task are lacking. To address this need, we propose a varying coefficient quantile regression model able to handle bivariate functional responses. Our work supports innovative insights into biomedical data by modeling the joint distribution of functional variables over their domains and across clinical covariates. We propose an estimation procedure based on the alternating direction method of multipliers and propagation separation algorithms to estimate varying coefficients using a B-spline basis and an $L_2$ smoothness penalty that encourages interpretability. A simulation study and an application to a real-world neurodevelopmental data set demonstrates the performance of our model and the insights provided by modeling functional fractional anisotropy and mean diffusivity jointly and their association with gestational age and sex.

A comprehensive comparison on cell-type composition inference for spatial transcriptomics data.

Spatial transcriptomics (ST) technologies allow researchers to examine transcriptional profiles along with maintained positional information. Such spatially resolved transcriptional characterization of intact tissue samples provides an integrated view of gene expression in its natural spatial and functional context. However, high-throughput sequencing-based ST technologies cannot yet reach single cell resolution. Thus, similar to bulk RNA-seq data, gene expression data at ST spot-level reflect transcriptional profiles of multiple cells and entail the inference of cell-type composition within each ST spot for valid and powerful subsequent analyses. Realizing the critical importance of cell-type decomposition, multiple groups have developed ST deconvolution methods. The aim of this work is to review state-of-the-art methods for ST deconvolution, comparing their strengths and weaknesses. In particular, we construct ST spots from single-cell level ST data to assess the performance of 10 methods, with either ideal reference or non-ideal reference. Furthermore, we examine the performance of these methods on spot- and bead-level ST data by comparing estimated cell-type proportions to carefully matched single-cell ST data. In comparing the performance on various tissues and technological platforms, we concluded that RCTD and stereoscope achieve more robust and accurate inferences.

Merging or ensembling: integrative analysis in multiple neuroimaging studies.

The aim of this paper is to systematically investigate merging and ensembling methods for spatially varying coefficient mixed effects models (SVCMEM) in order to carry out integrative learning of neuroimaging data obtained from multiple biomedical studies. The "merged" approach involves training a single learning model using a comprehensive dataset that encompasses information from all the studies. Conversely, the "ensemble" approach involves creating a weighted average of distinct learning models, each developed from an individual study. We systematically investigate the prediction accuracy of the merged and ensemble learners under the presence of different degrees of interstudy heterogeneity. Additionally, we establish asymptotic guidelines for making strategic decisions about when to employ either of these models in different scenarios, along with deriving optimal weights for the ensemble learner. To validate our theoretical results, we perform extensive simulation studies. The proposed methodology is also applied to 3 large-scale neuroimaging studies.

Urinary metabolite concentrations of phthalate and plasticizers in infancy and childhood in the UNC Baby Connectome Project.

INTRODUCTION: Existing evidence suggests that exposure to phthalates is higher among younger age groups. However, limited knowledge exists on how phthalate exposure, as well as exposure to replacement plasticizers, di(isononyl) cyclohexane-1,2-dicarboxylate (DINCH) and di-2-ethylhexyl terephthalate (DEHTP), change from infancy through early childhood. METHODS: Urine samples were collected across the first 5 years of life from typically developing infants and young children enrolled between 2017 and 2020 in the longitudinal UNC Baby Connectome Project. From 438 urine samples among 187 participants, we quantified concentrations of monobutyl phthalate (MnBP), mono-3-carboxypropyl phthalate (MCPP), monoisobutyl phthalate (MiBP), monoethyl phthalate (MEP) monobenzyl phthalate (MBzP), and metabolites of di(2-ethylhexyl) phthalate (DEHP), diisonoyl phthalate (DiNP), DINCH and DEHTP. Specific gravity (SG) adjusted metabolite and molar sum concentrations were compared across age groups. Intraclass correlation coefficients (ICCs) were calculated among 122 participants with multiple urine specimens (373 samples). RESULTS: Most phthalate metabolites showed high detection frequencies (>80% of samples). Replacement plasticizers DINCH (58-60%) and DEHTP (>97%) were also commonly found. DiNP metabolites were less frequently detected (<10%). For some metabolites, SG-adjusted concentrations were inversely associated with age, with the highest concentrations found in the first year of life. ICCs revealed low to moderate reliability in metabolite measurements (ρ = 0.10-0.48) suggesting a high degree of within-individual variation in exposure among this age group. The first 6 months (compared to remaining age groups) showed an increased ratio of carboxylated metabolites of DEHP and DEHTP, compared to other common metabolites, but no clear age trends for DINCH metabolite ratios. CONCLUSION: Metabolites of phthalates and replacements plasticizers were widely detected in infancy and early childhood, with the highest concentrations observed in the first year of life for several metabolites. Higher proportions of carboxylated metabolites of DEHP and DEHTP in younger age groups indicate potential differences in metabolism during infancy.

Addressing multi-site functional MRI heterogeneity through dual-expert collaborative learning for brain disease identification.

Several studies employ multi-site rs-fMRI data for major depressive disorder (MDD) identification, with a specific site as the to-be-analyzed target domain and other site(s) as the source domain. But they usually suffer from significant inter-site heterogeneity caused by the use of different scanners and/or scanning protocols and fail to build generalizable models that can well adapt to multiple target domains. In this article, we propose a dual-expert fMRI harmonization (DFH) framework for automated MDD diagnosis. Our DFH is designed to simultaneously exploit data from a single labeled source domain/site and two unlabeled target domains for mitigating data distribution differences across domains. Specifically, the DFH consists of a domain-generic student model and two domain-specific teacher/expert models that are jointly trained to perform knowledge distillation through a deep collaborative learning module. A student model with strong generalizability is finally derived, which can be well adapted to unseen target domains and analysis of other brain diseases. To the best of our knowledge, this is among the first attempts to investigate multi-target fMRI harmonization for MDD diagnosis. Comprehensive experiments on 836 subjects with rs-fMRI data from 3 different sites show the superiority of our method. The discriminative brain functional connectivities identified by our method could be regarded as potential biomarkers for fMRI-related MDD diagnosis.

Bayesian sparse heritability analysis with high-dimensional neuroimaging phenotypes.

Heritability analysis plays a central role in quantitative genetics to describe genetic contribution to human complex traits and prioritize downstream analyses under large-scale phenotypes. Existing works largely focus on modeling single phenotype and currently available multivariate phenotypic methods often suffer from scaling and interpretation. In this article, motivated by understanding how genetic underpinning impacts human brain variation, we develop an integrative Bayesian heritability analysis to jointly estimate heritabilities for high-dimensional neuroimaging traits. To induce sparsity and incorporate brain anatomical configuration, we impose hierarchical selection among both regional and local measurements based on brain structural network and voxel dependence. We also use a nonparametric Dirichlet process mixture model to realize grouping among single nucleotide polymorphism-associated phenotypic variations, providing biological plausibility. Through extensive simulations, we show the proposed method outperforms existing ones in heritability estimation and heritable traits selection under various scenarios. We finally apply the method to two large-scale imaging genetics datasets: the Alzheimer's Disease Neuroimaging Initiative and United Kingdom Biobank and show biologically meaningful results.

Altered neural flexibility in children with attention-deficit/hyperactivity disorder.

Attention-deficit/hyperactivity disorder (ADHD) is one of the most common neurodevelopmental disorders of childhood, and is often characterized by altered executive functioning. Executive function has been found to be supported by flexibility in dynamic brain reconfiguration. Thus, we applied multilayer community detection to resting-state fMRI data in 180 children with ADHD and 180 typically developing children (TDC) to identify alterations in dynamic brain reconfiguration in children with ADHD. We specifically evaluated MR derived neural flexibility, which is thought to underlie cognitive flexibility, or the ability to selectively switch between mental processes. Significantly decreased neural flexibility was observed in the ADHD group at both the whole brain (raw p = 0.0005) and sub-network levels (p < 0.05, FDR corrected), particularly for the default mode network, attention-related networks, executive function-related networks, and primary networks. Furthermore, the subjects with ADHD who received medication exhibited significantly increased neural flexibility (p = 0.025, FDR corrected) when compared to subjects with ADHD who were medication naïve, and their neural flexibility was not statistically different from the TDC group (p = 0.74, FDR corrected). Finally, regional neural flexibility was capable of differentiating ADHD from TDC (Accuracy: 77% for tenfold cross-validation, 74.46% for independent test) and of predicting ADHD severity using clinical measures of symptom severity (R2: 0.2794 for tenfold cross-validation, 0.156 for independent test). In conclusion, the present study found that neural flexibility is altered in children with ADHD and demonstrated the potential clinical utility of neural flexibility to identify children with ADHD, as well as to monitor treatment responses and disease severity.

Cross-trait prediction accuracy of summary statistics in genome-wide association studies.

In the era of big data, univariate models have widely been used as a workhorse tool for quickly producing marginal estimators; and this is true even when in a high-dimensional dense setting, in which many features are "true," but weak signals. Genome-wide association studies (GWAS) epitomize this type of setting. Although the GWAS marginal estimator is popular, it has long been criticized for ignoring the correlation structure of genetic variants (i.e., the linkage disequilibrium [LD] pattern). In this paper, we study the effects of LD pattern on the GWAS marginal estimator and investigate whether or not additionally accounting for the LD can improve the prediction accuracy of complex traits. We consider a general high-dimensional dense setting for GWAS and study a class of ridge-type estimators, including the popular marginal estimator and the best linear unbiased prediction (BLUP) estimator as two special cases. We show that the performance of GWAS marginal estimator depends on the LD pattern through the first three moments of its eigenvalue distribution. Furthermore, we uncover that the relative performance of GWAS marginal and BLUP estimators highly depends on the ratio of GWAS sample size over the number of genetic variants. Particularly, our finding reveals that the marginal estimator can easily become near-optimal within this class when the sample size is relatively small, even though it ignores the LD pattern. On the other hand, BLUP estimator has substantially better performance than the marginal estimator as the sample size increases toward the number of genetic variants, which is typically in millions. Therefore, adjusting for the LD (such as in the BLUP) is most needed when GWAS sample size is large. We illustrate the importance of our results by using the simulated data and real GWAS.

Asynchronous functional linear regression models for longitudinal data in reproducing kernel Hilbert space.

Motivated by the analysis of longitudinal neuroimaging studies, we study the longitudinal functional linear regression model under asynchronous data setting for modeling the association between clinical outcomes and functional (or imaging) covariates. In the asynchronous data setting, both covariates and responses may be measured at irregular and mismatched time points, posing methodological challenges to existing statistical methods. We develop a kernel weighted loss function with roughness penalty to obtain the functional estimator and derive its representer theorem. The rate of convergence, a Bahadur representation, and the asymptotic pointwise distribution of the functional estimator are obtained under the reproducing kernel Hilbert space framework. We propose a penalized likelihood ratio test to test the nullity of the functional coefficient, derive its asymptotic distribution under the null hypothesis, and investigate the separation rate under the alternative hypotheses. Simulation studies are conducted to examine the finite-sample performance of the proposed procedure. We apply the proposed methods to the analysis of multitype data obtained from the Alzheimer's Disease Neuroimaging Initiative (ADNI) study, which reveals significant association between 21 regional brain volume density curves and the cognitive function. Data used in preparation of this paper were obtained from the ADNI database (adni.loni.usc.edu).

Early brain development in infants at high risk for autism spectrum disorder.

Brain enlargement has been observed in children with autism spectrum disorder (ASD), but the timing of this phenomenon, and the relationship between ASD and the appearance of behavioural symptoms, are unknown. Retrospective head circumference and longitudinal brain volume studies of two-year olds followed up at four years of age have provided evidence that increased brain volume may emerge early in development. Studies of infants at high familial risk of autism can provide insight into the early development of autism and have shown that characteristic social deficits in ASD emerge during the latter part of the first and in the second year of life. These observations suggest that prospective brain-imaging studies of infants at high familial risk of ASD might identify early postnatal changes in brain volume that occur before an ASD diagnosis. In this prospective neuroimaging study of 106 infants at high familial risk of ASD and 42 low-risk infants, we show that hyperexpansion of the cortical surface area between 6 and 12 months of age precedes brain volume overgrowth observed between 12 and 24 months in 15 high-risk infants who were diagnosed with autism at 24 months. Brain volume overgrowth was linked to the emergence and severity of autistic social deficits. A deep-learning algorithm that primarily uses surface area information from magnetic resonance imaging of the brain of 6-12-month-old individuals predicted the diagnosis of autism in individual high-risk children at 24 months (with a positive predictive value of 81% and a sensitivity of 88%). These findings demonstrate that early brain changes occur during the period in which autistic behaviours are first emerging.

The emergence of a functionally flexible brain during early infancy.

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