Neuroanatomical, transcriptomic, and molecular correlates of math ability and their prognostic value for predicting learning outcomes
In our increasingly technology-driven world, mathematical thinking underpins virtually every aspect of our daily lives. Foundational mathematical abilities, acquired in early childhood, are prerequisites for individual academic and professional development and success . Despite their ubiquity and importance, the neurobiological mechanisms that drive some children to have better (or poorer) mathematical abilities and learning outcomes than others have rarely been investigated. A comprehensive understanding of the macroscale structure of the brain and the microscale transcriptomic/molecular factors and their interplay, particularly in relation to mathematical abilities and learning, is currently lacking. This knowledge is critical for the comprehensive multimodal characterization of the neurobiological mechanisms driving individual differences in mathematical abilities and learning outcomes. To address this knowledge gap, our study aims to integrate neuroimaging and transcriptomic data to probe the neurobiological mechanisms underlying mathematical abilities and learning in children. Neuroimaging offers insights into the structural and functional brain patterns associated with mathematical cognition, while genetic analysis reveals the underlying transcriptomic and molecular factors. By combining these methodologies and applying them to multiple datasets ( and figs. S1 and S2), we sought to elucidate the interrelations between brain structure, cognitive function, genetics, and physiology and enhance our understanding of individual differences in mathematical cognition and learning. This integrated approach is not only crucial for advancing scientific knowledge but also has substantial implications for precision education. We focused our investigation on three targeted questions in school-age children: (i) What are the replicable brain structural correlates of individual differences in mathematical abilities? (ii) Which transcriptomic molecular factors contribute to this brain-behavior association? (iii) Can transcriptomic molecular factors predict learning outcomes in response to math interventions? By addressing these questions, we aim to lay the groundwork for a comprehensive multilevel characterization of brain and transcriptomic and molecular correlates of children’s mathematical abilities and learning during a pivotal stage of knowledge acquisition. The first aim of our study was to identify brain structural correlates of individual differences in children’s mathematical abilities. A substantial body of research has shown that mathematical cognition relies on multiple distributed brain systems, including the bilateral posterior parietal cortex (PPC), ventral temporal-occipital cortex (VTOC), and prefrontal cortex (PFC). Collectively, these brain regions support quantity representation, symbolic number form, as well as higher-order functions such as working memory and cognitive control required for efficient problem-solving . While prior studies have examined the neuroanatomical sources of differences in children’s mathematical abilities, no consensus has yet emerged (table S1) . Crucially, the reliance on categorical analyses comparing atypically and typically developing children in many prior studies, along with the use of arbitrary cutoffs on ability measures, has made it difficult to draw definitive conclusions regarding the relationship between brain structure and individual differences in children’s mathematical abilities. Findings from a few studies with dimensional approaches that investigated the relation between brain structure and individual differences in mathematical abilities in children have also been largely inconsistent . Furthermore, to our knowledge, the replicability of these findings across multiple cohorts remains unknown. The lack of consistent findings regarding brain structural correlates of individual differences in children’s mathematical abilities likely stems from small sample sizes in previous research (see table S1 for a summary). Specifically, the small sample sizes ( n < 25 in each group) in many of these studies have hindered the ability to detect robust associations with gray matter volume and morphology. Moreover, most investigations have relied on univariate brain-behavior associations, linking a single brain region to a mathematical ability measure. Such a regional approach is inconsistent with the systems neuroscience perspective that multiple distributed brain areas contribute to complex cognitive processes involved in mathematical reasoning . In addition, univariate approaches, compared to multivariate methods, have been shown to be less reliable in predicting brain-behavior relationships . We used a multivariate analysis strategy, leveraging two large cohorts of structural magnetic resonance imaging (MRI) data combined with standardized neuropsychological assessments of mathematical ability . Our study included a primary discovery cohort comprising data acquired at Stanford University, consisting of 219 children of ages between 7 and 13 years, and an independent Nathan Kline Institute-Rockland Sample (NKI-RS) validation cohort, which included data from 91 children within the same age range ( and fig. S1). We used canonical correlation analysis, a widely used multivariate analysis approach in brain-behavior association studies, to uncover associations between brain and behavioral measures by identifying maximally correlated linear combinations between two sets of variables . Our study features one of the largest cohorts of brain structure–mathematical abilities data, specifically targeting school-age children during a critical period of knowledge acquisition. To our knowledge, this represents a unique multicohort dataset of its kind in the field. We hypothesized that multivariate patterns of structural brain features, encompassing the posterior parietal, ventral temporal-occipital, and prefrontal cortical areas, would predict individual differences in children’s mathematical abilities. We further hypothesized that these associations would be replicated across the two independent cohorts, thereby establishing a mathematical ability–related imaging phenotype (MAIP). The second aim of our study was to identify transcriptomic molecular factors associated with the MAIP. Genetic factors are thought to play a critical role in the development of brain structures that support complex cognitive functions such as mathematics . Prior research has suggested that mathematical abilities may be highly heritable, with genetic factors accounting for 60 to 70% of the observed variance in twin studies . One recent study pinpointed a specific gene, ROBO1 , that is significantly associated with gray matter volume in the right inferior parietal cortex, a region implicated in mathematical cognition . However, it is unlikely that complex cognitive functions like mathematical abilities are solely supported by single gene–brain region association . Notably, a groundbreaking genome-wide association study (GWAS) by Lee et al.  demonstrated that the genetic landscape is multifaceted, identifying 618 and 365 single-nucleotide polymorphisms (SNPs) associated with mathematical ability and the highest math class taken, respectively, from a sample of 1.1 million individuals. We examined how genome-wide gene expression, encompassing known candidate genes, varies across brain regions linked to mathematical abilities. Our approach offers a broader and more detailed perspective on transcriptomic and molecular factors linked to mathematical cognition, thereby enhancing our comprehension of the biological underpinnings of these abilities. The Allen Human Brain Atlas, a comprehensive genome-wide atlas of gene expression built on ~3000 tissue samples covering the whole human brain with high spatial resolution , provides an unprecedented opportunity to bridge the gap between brain imaging and transcriptome measures . This atlas addresses the limitations of previous studies , which lack precise information regarding variations in gene-expression profiles across individual brain regions . In recent years, the Allen Human Brain Atlas has been instrumental in elucidating the transcriptomic and molecular factors that shape large-scale human brain organization, both in neurotypical adults , as well as normative and atypical brain development . These studies have unveiled a close link between the spatial patterns of gene expression and brain structure, indicating that transcriptional signatures can capture the underlying variations in molecular, biological, and cellular functions that drive functional specialization in the human brain . We leveraged the Allen Human Brain Atlas and conducted partial least squares analysis to probe connections between gene expression and brain structures that account for individual differences in children’s mathematical abilities . Our analysis yielded two key measures: (i) brain-wide gene expression patterns related to spatial variability in MAIP and (ii) an annotated list of genes ranked by their impact on MAIP. We hypothesized that specific patterns of brain-wide gene expression would explain MAIP and that the most influential genes in these expression profiles would overlap with previously reported candidate genes related to mathematical abilities . Furthermore, we used BrainSpan to examine the developmental trajectory of the observed gene expression profiles associated with MAIP. The third aim of our study was to investigate whether transcriptomic molecular factors correlated with the MAIP could serve as reliable predictive markers for individual differences in mathematical learning outcomes following cognitive interventions. Previous research has shown that while early math interventions are effective for many children who experience difficulties with mathematics , learning outcomes following the interventions have been shown to vary considerably across individuals . A particularly intriguing aspect of academic performance is the relatively high stability of mathematical abilities in children throughout their school years. Critically, the stability of educational and mathematical achievement across school years is shown to be largely explained by genetic factors, accounting for about 60% of the stability in educational achievement from early childhood to late adolescence . This observation not only emphasizes the essential role of genetics but also underscores the critical need to understand the underlying mechanisms that shape the interplay between genes and educational outcomes. In the present study, we sought to address this gap by examining the prognostic value of transcriptomic molecular factors associated with the MAIP for predicting learning outcomes following cognitive interventions designed to improve mathematical skills. We leveraged brain imaging and learning measures from two distinct math interventions, in which children received training designed to enhance one of two foundational numerical skills: (i) retrieval of math facts  or (ii) number sense . The two types of interventions were designed to remediate poor mathematical abilities in children, including those with mathematical learning disabilities. Children in these tutoring programs have demonstrated considerable variability in math learning outcomes in response to the interventions . Transcriptome-imaging studies have begun to open avenues for the development of individual-level predictive tools. A notable advance in this field was made by Di Biase and colleagues , who developed an individual-level transcriptome similarity index (TSI) by correlating neuroimaging-based brain deviation maps of patients with the gene expression maps of specific cell classes. This work demonstrated the functional relevance of the TSI by linking it with individual polygenic cell-class scores derived from genotype data in patients with schizophrenia. Building on this innovative and validated approach, we used similar methodologies to identify transcriptome-based predictors for individual learning variations in response to math interventions. We developed a TSI to predict math learning outcomes; the TSI was derived by correlating the transcriptomic molecular factors associated with the MAIP with the individual brain structure patterns of each child. This index was calculated for every participant in the two math intervention programs, using structural MRI data collected before the commencement of tutoring. We hypothesized that TSI would consistently predict mathematical learning outcomes across diverse learning contexts, thereby offering important insights into the fundamental role of gene expression in learning and skill acquisition. Our study illuminates the complex interplay between brain structure, gene expression, and educational achievement and uncovers transcriptome-based predictors to enable the identification of children who may benefit from more intensive intervention. Our findings have the potential to substantially enhance the development of personalized educational approaches.   Identifying MAIP: Association between individual differences in mathematical ability and brain structure We first used the Stanford discovery cohort to determine the relationship between individual differences in mathematical ability and brain structure. Behavioral and structural MRI brain imaging data from 219 children (age: 7 to 13 years, M = 9.8 years, SD = 1.7 years; 99 females; ) who met criteria for high-quality data were used in the analysis (see Materials and Methods). Mathematical abilities were measured using grade-normed standardized scores of numerical operations and math reasoning subtests from the Wechsler Individual Achievement Test Second Edition (WIAT-II) . The numerical operations subtest assessed number writing, identification, production, and simple arithmetic skills, while the math reasoning subtest assessed counting, geometric shape identification, and math problem-solving (numerical operations: M = 110.5, SD = 20.1; math reasoning: M = 113.2, SD = 15.6; ). Gray matter volume was computed in each of the 246 regions of interest defined by the Brainnetome atlas . Subsequently, a multivariate canonical correlation analysis  was conducted across 219 children, with brain measures comprising gray matter volumes of the 246 regions of interest  and behavioral measures including WIAT-II numerical operations and math reasoning scores. The canonical correlation analysis revealed a significant relationship between gray matter volume and mathematical abilities ( r = 0.34, P = 0.023, permutation test 1000 times; ). illustrates the canonical correlation weights for brain and behavioral measures. Notably, the brain pattern/brain measure weight identified in the canonical correlation analysis was correlated with both math measures, WIAT-II numerical operations and math reasoning (with correlation coefficients of −0.7088 and −0.9999, respectively), likely reflecting a shared neuroanatomical basis for mathematical abilities across the two assessments. The results demonstrate that individual differences in mathematical abilities in children are associated with patterns of widespread variations in gray matter volume. Specifically, poorer (better) performance on both WIAT-II numerical operations and math reasoning was associated with higher (lower) gray matter volume in PPC, VTOC, and PFC—regions implicated in numerical cognition —and lower (higher) gray matter volume in other regions, including visual, posterior insular, and subcortical regions. Additional control analyses revealed a similar canonical correlation mode under different analysis strategies including permutation inference , controlling for potential confounds, such as age and sex, and dimension reduction thresholds (Supplementary Text and table S2), highlighting the robustness of our findings. Furthermore, we validated our findings by implementing a 10-fold cross-validation approach and testing prediction of mathematical ability from an alternative assessment. Significant results were obtained from both the analyses, affirming the generalizability of our findings across participants and math assessments (Supplementary Text and fig. S3). In addition, no significant prediction was observed for nonmath measures including working memory or socioeconomic status (fig. S4), suggesting that the observed brain structure patterns were specifically associated with mathematical ability rather than broadly representing general cognitive abilities or other noncognitive factors. Together, these results demonstrate that individual differences in mathematical abilities are associated with distinct patterns of higher and lower gray matter volume distributed across the brain and provide an index for a MAIP. Association between MAIP and math-related meta-analytic brain activation maps To confirm the cognitive relevance of the identified MAIP pattern, we compared it with math-related brain functional activation maps available in the Neurosynth database derived from term-based meta-analyses . The math-related terms in the Neurosynth database included “symbolic,” “arithmetic,” “subtraction,” “calculation,” and “addition,” pinpointing multiple brain regions consistently activated during math task performance in functional neuroimaging studies . We performed a spatial correlation analysis between the identified MAIP pattern and the meta-analytic brain activation map for each term. We found that the regional brain activation distribution in all of these math-related terms was significantly correlated to the MAIP pattern ( P s < 0.001, permutation test 1000 times; ). Notably, the brain regions identified in these math-related terms included areas with higher gray matter volume in the MAIP, including the PPC, VTOC, and PFC. Similar spatial correlation ( P perm < 0.001; fig. S5) was observed between the identified MAIP pattern and a math-related brain activation map derived from a meta-analysis of math-related functional neuroimaging studies exclusively focused on children (table S3). Furthermore, we observed that the spatial correlations between the MAIP and meta-analytic brain activation maps for math-related terms were significantly higher than those between the MAIP and reading-related control terms, such as reading and writing ( P s perm < 0.001; and table S4), implying cognitive specificity for the spatial distribution of the MAIP pattern. Similar results were also observed in additional analyses using “math” and “reading” topic maps instead of term maps in the NeuroSynth database (fig. S6). These results demonstrate that the identified MAIP effectively captures brain organization and regional specialization specifically related to mathematical information processing, corroborating findings from task-based functional neuroimaging studies. Generalizing MAIP in an independent cohort: Replication of association between individual differences in mathematical ability and brain structure Next, we sought to investigate the generalizability of the findings drawn from the Stanford discovery cohort by examining an independent cohort whose age range matched that of the Stanford cohort. Specifically, we used data from the NKI-RS, a publicly available life-span dataset consisting of a large-scale community sample in Rockland County . Ninety-one children with high-quality structural brain MRI data and numerical operations scores from the WIAT-II were included in the NKI-RS cohort (age: M = 10.1 years, SD = 1.9 years; 45 females; and ). Math reasoning scores from the WIAT-II were not available in the NKI-RS cohort. We computed gray matter volumes of 246 regions of interest for each child . To predict WIAT-II numerical operations scores for children within the NKI-RS cohort, we developed a prediction approach that involved multiplying the gray matter volume of each child by the canonical vector of the observed canonical correlation mode from the Stanford cohort. This process generated a math-weighted brain score, which was then adjusted by scaling it using the mean and SD of WIAT-II numerical operation scores from the Stanford cohort (see Materials and Methods for details). The correlation between the predicted and observed WIAT-II numerical operations scores across the children in the NKI-RS cohort was then calculated to evaluate the accuracy of prediction. We found a significant correlation between the predicted and observed WIAT-II numerical operations scores within this independent validation cohort ( r = 0.29, P = 0.005; ). Notably, no significant correlation was observed between the predicted score and word reading scores from the WIAT-II within the NKI-RS cohort ( r = 0.07, P = 0.484). Furthermore, a direct comparison of correlation coefficients for math and reading ability revealed a significant difference between these two domains ( t = 2.28, P = 0.025, Williams’s test), suggesting that the generalizability of our findings from the Stanford cohort pertains specifically to mathematical ability. Additional control analyses with other nonmathematical measures reinforced the specificity of math ability–related findings (Supplementary Text). These results provide evidence of a generalizable MAIP pattern for specifically predicting mathematical ability in an independent cohort of age-matched children. Brain gene expression profiles associated with the MAIP We examined gene expression profiles associated with the MAIP pattern using human brain-wide gene expression data from the Allen Human Brain Atlas . Gene expression maps, consisting of 246 regions and 15,633 genes, were created through preprocessing using standardized imaging-transcriptomics methodology . We investigated the relationship between brain-wide gene expression maps and the MAIP using a partial least squares regression analysis . The partial least squares regression analysis revealed that the top three significant components explained 72% of the variance in the MAIP pattern ( P s < 0.001, permutation test 1000 times, corrected for spatial autocorrelation; and fig. S7) . Notably, the first component explained more than 50% of its spatial variance, while the second and third components each explained around 10% of the variance, pointing to the dominant role of the first component. This dominant first component was positively correlated with the MAIP and was characterized by a high gene expression profile in the parietal-prefrontal areas, including PPC, PFC, as well as the right VTOC , consistent with brain areas that showed higher gray matter volume in the MAIP. The other two components identified from the partial least squares regression are described in Supplementary Text and fig. S8. In addition, the partial least squares regression analysis yielded an ordered gene list for each significant component, based on the weight of genes contributing to the gene expression profiles. Information about three ordered gene lists is available in Data and materials availability. Together, these results suggest that distinct brain-wide gene expression profiles underlie the MAIP pattern, corresponding to the distributed pattern of varying gray matter volume underlying individual differences in mathematical abilities. Developmental trajectory of brain gene expression profiles associated with the MAIP We subsequently examined the developmental trajectory of gene expression profiles associated with the MAIP using human developmental gene expression data from BrainSpan . Gene expression data from 469 brain tissue samples spanning 15 brain regions from participants aged from eight postconceptional weeks to 40 years old were downloaded. Genes that overlapped with the Allen Human Brain Atlas analyses were examined. Brain tissue samples were categorized into six developmental periods: prenatal, infancy, early childhood, late childhood, adolescence, and adulthood. Notably, the age range of the late childhood group overlapped with the age range of our sample. We then calculated the MAIP-weighted average gene expression within each age bin for each brain region made available through BrainSpan. This yielded a developmental trajectory of brain region-wise gene expression profiles associated with the MAIP. Briefly, the developmental trajectory is indicative of how expression profiles of genes associated with mathematical ability in a brain region change with age. Notably, as shown in fig. S9, the gene expression related to the MAIP remained relatively high and stable from late childhood to adulthood across multiple brain regions. These results indicate the enduring stability of MAIP-related gene expression profiles throughout later stages of development, which aligns with earlier observations that underscore of educational achievement stability persisting from childhood to adolescence . Association of MAIP-related gene lists with candidate math ability–related genes using gene set enrichment analysis Next, we conducted gene set enrichment analysis to determine whether previously GWAS-reported candidate genes related to math are overrepresented in the most strongly correlated genes identified in the ordered gene lists associated with the MAIP. The candidate math ability–related genes were selected on the basis of the largest GWAS to date, involving 1.1 million individuals , which identified SNPs located in 386 genes associated with mathematical ability and 248 genes associated with the highest math class taken. We used the ordered gene lists in descending order for each of the three significant components derived from the partial least squares analysis. For the dominant first component, we found a significant enrichment of GWAS-identified gene sets related to mathematical ability and the highest math class taken at the top of the ordered gene list. In other words, the genes most positively correlated with the weighted gene expression profile associated with the MAIP were overrepresented by candidate math ability–related genes [math ability: normalized enrichment score ( NES ) = 1.85, P adjust < 0.001; highest math class taken: NES = 1.62, P adjust = 0.002; and table S5]. Additional partial least squares regression analysis revealed a similar gene expression profile pattern for MAIP when considering only the math-related genes (Supplementary Text and fig. S10). These results confirm the critical role of the candidate math ability–related genes in the expression profiles associated with the MAIP pattern. We additionally examined whether previously GWAS-reported candidate genes related to reading-related abilities and disabilities, i.e., dyslexia , exhibited any association with the MAIP. We found that gene sets related to reading did not show significant enrichment within our ordered gene list of the first dominant component [ NES = −0.89, P adjust > 0.898 for the gene list from Doust et al. , ; NES = 0.64, P adjust > 0.898 for the gene list from ; NES = −0.68, P adjust > 0.898 for the gene list from Price et al. ; fig. S11]. Additional analysis showed that genes most positively correlated with the MAIP were not significantly overrepresented by working memory–related genes  ( NES = 1.54, P adjust = 0.100; fig. S12), which suggests distinct profiles for math ability–related genes and genes related to more general cognitive abilities. Findings of the other two components are described in Supplementary Text. These results reveal that MAIP-related genes are specifically overrepresented by candidate math ability–related genes, but not reading or working memory–related genes, implying a distinct link between MAIP and the up-regulation of candidate math ability–related genes. Functional characterization of MAIP-related gene lists using Gene Ontology enrichment analysis To elucidate the functional characteristics of the MAIP genes, we conducted a Gene Ontology enrichment analysis. All three Gene Ontology classes were included: biological process, cellular component, and molecular function. For the dominant first component, we found that the top of the descending-ordered gene list, i.e., the genes most positively related to the weighted gene expression profile, were significantly enriched for genes related to neurons, neuronal signaling and synaptic transmission, and voltage-gated potassium channel activity [false discovery rate (FDR)–corrected q s < 0.001; and figs. S13 to S15; see figs. S16 and S17 for the genes most negatively related to this profile]. Functional characteristics of gene expression profiles for the other two components are described in Supplementary Text and figs. S18 to S21. These results reveal that the MAIP is associated with the transcriptional profiles of genes that are enriched for neurons and voltage-gated potassium channel activity. Linking MAIP-related gene expression profiles and brain structure to children’s math learning outcomes To determine the prognostic value of the observed MAIP-related gene expression profile, we developed a TSI for the MAIP gene expression profile with each child’s brain structure and determined whether this similarity index can predict children’s math learning outcomes. We first examined an 8-week intervention dataset comprising 24 children (age: 7 to 10 years, M = 8.4 years, SD = 0.5 years; 14 females; ) . Gray matter volume of 246 brain regions was estimated for each child on the basis of the structural brain image acquired before tutoring. Learning outcome was measured as the improvement in processing speed for math problems in an arithmetic verification task in response to tutoring (Materials and Methods). We computed three TSIs in each child, corresponding to each of the three significant gene expression profiles identified by the partial least squares regression analysis. Each index was calculated as the correlation between MAIP gene expression profiles and gray matter volumes in each child. We then related these indices to learning outcomes and generated predicted scores using a general linear model, with TSIs as predictor variables and learning outcome as the response variable. As shown in , these three similarity indices for the three components together significantly predicted learning outcomes in the 8-week intervention dataset ( r = 0.62, P < 0.001). Similar results were observed when using the similarity index of only the dominant first component, with a large TSI associated with intervention-induced improvement in response time for math problems ( r = −0.43, P = 0.034). These results remained significant when controlling for age or sex (Supplementary Text). We additionally examined whether behavioral or brain anatomical measures alone could predict children’s future learning outcomes in response to math interventions and found that our TSI, which combined multimodal information, outperformed unimodal behavioral and brain anatomical measures in predicting intervention-induced math learning (Supplementary Text). Linking MAIP-related gene expression profile and brain structure to children’s math learning outcomes: Replication and generalizability analysis Next, we sought to replicate and demonstrate the generalizability of these findings using a 4-week intervention dataset ( n = 61 children; age: M = 8.2 years, SD = 0.6 years; 33 females; ), which constituted an independent intervention sample that also did not overlap with the cohort used to construct the MAIP . Gray matter volume of 246 brain regions was estimated for each child on the basis of the structural brain image acquired before tutoring. Learning outcome was measured as improvement in the numerical distance effect of processing speed for math problems in a symbolic quantity discrimination task in response to tutoring (Materials and Methods). As described previously, we computed three TSIs in each child and generated predicted scores using a general linear model. As shown in , these three similarity indices for the three components together significantly predicted learning outcomes in the 4-week intervention dataset ( r = 0.38, P = 0.003). Similar results were observed using the similarity index of only the dominant first component, with a large TSI associated with intervention-induced improvement in response time for math problems ( r = −0.25, P = 0.048). These results remained consistent when controlling for age or sex (Supplementary Text). Conventional unimodal prediction analysis yielded similar results as the 8-week intervention cohort (Supplementary Text). These results suggest that our multimodal TSI can reliably predict children’s learning outcomes in both intervention cohorts despite differences in the tutoring protocol, markedly outperforming conventional unimodal approaches.   We first used the Stanford discovery cohort to determine the relationship between individual differences in mathematical ability and brain structure. Behavioral and structural MRI brain imaging data from 219 children (age: 7 to 13 years, M = 9.8 years, SD = 1.7 years; 99 females; ) who met criteria for high-quality data were used in the analysis (see Materials and Methods). Mathematical abilities were measured using grade-normed standardized scores of numerical operations and math reasoning subtests from the Wechsler Individual Achievement Test Second Edition (WIAT-II) . The numerical operations subtest assessed number writing, identification, production, and simple arithmetic skills, while the math reasoning subtest assessed counting, geometric shape identification, and math problem-solving (numerical operations: M = 110.5, SD = 20.1; math reasoning: M = 113.2, SD = 15.6; ). Gray matter volume was computed in each of the 246 regions of interest defined by the Brainnetome atlas . Subsequently, a multivariate canonical correlation analysis  was conducted across 219 children, with brain measures comprising gray matter volumes of the 246 regions of interest  and behavioral measures including WIAT-II numerical operations and math reasoning scores. The canonical correlation analysis revealed a significant relationship between gray matter volume and mathematical abilities ( r = 0.34, P = 0.023, permutation test 1000 times; ). illustrates the canonical correlation weights for brain and behavioral measures. Notably, the brain pattern/brain measure weight identified in the canonical correlation analysis was correlated with both math measures, WIAT-II numerical operations and math reasoning (with correlation coefficients of −0.7088 and −0.9999, respectively), likely reflecting a shared neuroanatomical basis for mathematical abilities across the two assessments. The results demonstrate that individual differences in mathematical abilities in children are associated with patterns of widespread variations in gray matter volume. Specifically, poorer (better) performance on both WIAT-II numerical operations and math reasoning was associated with higher (lower) gray matter volume in PPC, VTOC, and PFC—regions implicated in numerical cognition —and lower (higher) gray matter volume in other regions, including visual, posterior insular, and subcortical regions. Additional control analyses revealed a similar canonical correlation mode under different analysis strategies including permutation inference , controlling for potential confounds, such as age and sex, and dimension reduction thresholds (Supplementary Text and table S2), highlighting the robustness of our findings. Furthermore, we validated our findings by implementing a 10-fold cross-validation approach and testing prediction of mathematical ability from an alternative assessment. Significant results were obtained from both the analyses, affirming the generalizability of our findings across participants and math assessments (Supplementary Text and fig. S3). In addition, no significant prediction was observed for nonmath measures including working memory or socioeconomic status (fig. S4), suggesting that the observed brain structure patterns were specifically associated with mathematical ability rather than broadly representing general cognitive abilities or other noncognitive factors. Together, these results demonstrate that individual differences in mathematical abilities are associated with distinct patterns of higher and lower gray matter volume distributed across the brain and provide an index for a MAIP.   To confirm the cognitive relevance of the identified MAIP pattern, we compared it with math-related brain functional activation maps available in the Neurosynth database derived from term-based meta-analyses . The math-related terms in the Neurosynth database included “symbolic,” “arithmetic,” “subtraction,” “calculation,” and “addition,” pinpointing multiple brain regions consistently activated during math task performance in functional neuroimaging studies . We performed a spatial correlation analysis between the identified MAIP pattern and the meta-analytic brain activation map for each term. We found that the regional brain activation distribution in all of these math-related terms was significantly correlated to the MAIP pattern ( P s < 0.001, permutation test 1000 times; ). Notably, the brain regions identified in these math-related terms included areas with higher gray matter volume in the MAIP, including the PPC, VTOC, and PFC. Similar spatial correlation ( P perm < 0.001; fig. S5) was observed between the identified MAIP pattern and a math-related brain activation map derived from a meta-analysis of math-related functional neuroimaging studies exclusively focused on children (table S3). Furthermore, we observed that the spatial correlations between the MAIP and meta-analytic brain activation maps for math-related terms were significantly higher than those between the MAIP and reading-related control terms, such as reading and writing ( P s perm < 0.001; and table S4), implying cognitive specificity for the spatial distribution of the MAIP pattern. Similar results were also observed in additional analyses using “math” and “reading” topic maps instead of term maps in the NeuroSynth database (fig. S6). These results demonstrate that the identified MAIP effectively captures brain organization and regional specialization specifically related to mathematical information processing, corroborating findings from task-based functional neuroimaging studies.   Next, we sought to investigate the generalizability of the findings drawn from the Stanford discovery cohort by examining an independent cohort whose age range matched that of the Stanford cohort. Specifically, we used data from the NKI-RS, a publicly available life-span dataset consisting of a large-scale community sample in Rockland County . Ninety-one children with high-quality structural brain MRI data and numerical operations scores from the WIAT-II were included in the NKI-RS cohort (age: M = 10.1 years, SD = 1.9 years; 45 females; and ). Math reasoning scores from the WIAT-II were not available in the NKI-RS cohort. We computed gray matter volumes of 246 regions of interest for each child . To predict WIAT-II numerical operations scores for children within the NKI-RS cohort, we developed a prediction approach that involved multiplying the gray matter volume of each child by the canonical vector of the observed canonical correlation mode from the Stanford cohort. This process generated a math-weighted brain score, which was then adjusted by scaling it using the mean and SD of WIAT-II numerical operation scores from the Stanford cohort (see Materials and Methods for details). The correlation between the predicted and observed WIAT-II numerical operations scores across the children in the NKI-RS cohort was then calculated to evaluate the accuracy of prediction. We found a significant correlation between the predicted and observed WIAT-II numerical operations scores within this independent validation cohort ( r = 0.29, P = 0.005; ). Notably, no significant correlation was observed between the predicted score and word reading scores from the WIAT-II within the NKI-RS cohort ( r = 0.07, P = 0.484). Furthermore, a direct comparison of correlation coefficients for math and reading ability revealed a significant difference between these two domains ( t = 2.28, P = 0.025, Williams’s test), suggesting that the generalizability of our findings from the Stanford cohort pertains specifically to mathematical ability. Additional control analyses with other nonmathematical measures reinforced the specificity of math ability–related findings (Supplementary Text). These results provide evidence of a generalizable MAIP pattern for specifically predicting mathematical ability in an independent cohort of age-matched children.   We examined gene expression profiles associated with the MAIP pattern using human brain-wide gene expression data from the Allen Human Brain Atlas . Gene expression maps, consisting of 246 regions and 15,633 genes, were created through preprocessing using standardized imaging-transcriptomics methodology . We investigated the relationship between brain-wide gene expression maps and the MAIP using a partial least squares regression analysis . The partial least squares regression analysis revealed that the top three significant components explained 72% of the variance in the MAIP pattern ( P s < 0.001, permutation test 1000 times, corrected for spatial autocorrelation; and fig. S7) . Notably, the first component explained more than 50% of its spatial variance, while the second and third components each explained around 10% of the variance, pointing to the dominant role of the first component. This dominant first component was positively correlated with the MAIP and was characterized by a high gene expression profile in the parietal-prefrontal areas, including PPC, PFC, as well as the right VTOC , consistent with brain areas that showed higher gray matter volume in the MAIP. The other two components identified from the partial least squares regression are described in Supplementary Text and fig. S8. In addition, the partial least squares regression analysis yielded an ordered gene list for each significant component, based on the weight of genes contributing to the gene expression profiles. Information about three ordered gene lists is available in Data and materials availability. Together, these results suggest that distinct brain-wide gene expression profiles underlie the MAIP pattern, corresponding to the distributed pattern of varying gray matter volume underlying individual differences in mathematical abilities.   We subsequently examined the developmental trajectory of gene expression profiles associated with the MAIP using human developmental gene expression data from BrainSpan . Gene expression data from 469 brain tissue samples spanning 15 brain regions from participants aged from eight postconceptional weeks to 40 years old were downloaded. Genes that overlapped with the Allen Human Brain Atlas analyses were examined. Brain tissue samples were categorized into six developmental periods: prenatal, infancy, early childhood, late childhood, adolescence, and adulthood. Notably, the age range of the late childhood group overlapped with the age range of our sample. We then calculated the MAIP-weighted average gene expression within each age bin for each brain region made available through BrainSpan. This yielded a developmental trajectory of brain region-wise gene expression profiles associated with the MAIP. Briefly, the developmental trajectory is indicative of how expression profiles of genes associated with mathematical ability in a brain region change with age. Notably, as shown in fig. S9, the gene expression related to the MAIP remained relatively high and stable from late childhood to adulthood across multiple brain regions. These results indicate the enduring stability of MAIP-related gene expression profiles throughout later stages of development, which aligns with earlier observations that underscore of educational achievement stability persisting from childhood to adolescence .   Next, we conducted gene set enrichment analysis to determine whether previously GWAS-reported candidate genes related to math are overrepresented in the most strongly correlated genes identified in the ordered gene lists associated with the MAIP. The candidate math ability–related genes were selected on the basis of the largest GWAS to date, involving 1.1 million individuals , which identified SNPs located in 386 genes associated with mathematical ability and 248 genes associated with the highest math class taken. We used the ordered gene lists in descending order for each of the three significant components derived from the partial least squares analysis. For the dominant first component, we found a significant enrichment of GWAS-identified gene sets related to mathematical ability and the highest math class taken at the top of the ordered gene list. In other words, the genes most positively correlated with the weighted gene expression profile associated with the MAIP were overrepresented by candidate math ability–related genes [math ability: normalized enrichment score ( NES ) = 1.85, P adjust < 0.001; highest math class taken: NES = 1.62, P adjust = 0.002; and table S5]. Additional partial least squares regression analysis revealed a similar gene expression profile pattern for MAIP when considering only the math-related genes (Supplementary Text and fig. S10). These results confirm the critical role of the candidate math ability–related genes in the expression profiles associated with the MAIP pattern. We additionally examined whether previously GWAS-reported candidate genes related to reading-related abilities and disabilities, i.e., dyslexia , exhibited any association with the MAIP. We found that gene sets related to reading did not show significant enrichment within our ordered gene list of the first dominant component [ NES = −0.89, P adjust > 0.898 for the gene list from Doust et al. , ; NES = 0.64, P adjust > 0.898 for the gene list from ; NES = −0.68, P adjust > 0.898 for the gene list from Price et al. ; fig. S11]. Additional analysis showed that genes most positively correlated with the MAIP were not significantly overrepresented by working memory–related genes  ( NES = 1.54, P adjust = 0.100; fig. S12), which suggests distinct profiles for math ability–related genes and genes related to more general cognitive abilities. Findings of the other two components are described in Supplementary Text. These results reveal that MAIP-related genes are specifically overrepresented by candidate math ability–related genes, but not reading or working memory–related genes, implying a distinct link between MAIP and the up-regulation of candidate math ability–related genes.   To elucidate the functional characteristics of the MAIP genes, we conducted a Gene Ontology enrichment analysis. All three Gene Ontology classes were included: biological process, cellular component, and molecular function. For the dominant first component, we found that the top of the descending-ordered gene list, i.e., the genes most positively related to the weighted gene expression profile, were significantly enriched for genes related to neurons, neuronal signaling and synaptic transmission, and voltage-gated potassium channel activity [false discovery rate (FDR)–corrected q s < 0.001; and figs. S13 to S15; see figs. S16 and S17 for the genes most negatively related to this profile]. Functional characteristics of gene expression profiles for the other two components are described in Supplementary Text and figs. S18 to S21. These results reveal that the MAIP is associated with the transcriptional profiles of genes that are enriched for neurons and voltage-gated potassium channel activity.   To determine the prognostic value of the observed MAIP-related gene expression profile, we developed a TSI for the MAIP gene expression profile with each child’s brain structure and determined whether this similarity index can predict children’s math learning outcomes. We first examined an 8-week intervention dataset comprising 24 children (age: 7 to 10 years, M = 8.4 years, SD = 0.5 years; 14 females; ) . Gray matter volume of 246 brain regions was estimated for each child on the basis of the structural brain image acquired before tutoring. Learning outcome was measured as the improvement in processing speed for math problems in an arithmetic verification task in response to tutoring (Materials and Methods). We computed three TSIs in each child, corresponding to each of the three significant gene expression profiles identified by the partial least squares regression analysis. Each index was calculated as the correlation between MAIP gene expression profiles and gray matter volumes in each child. We then related these indices to learning outcomes and generated predicted scores using a general linear model, with TSIs as predictor variables and learning outcome as the response variable. As shown in , these three similarity indices for the three components together significantly predicted learning outcomes in the 8-week intervention dataset ( r = 0.62, P < 0.001). Similar results were observed when using the similarity index of only the dominant first component, with a large TSI associated with intervention-induced improvement in response time for math problems ( r = −0.43, P = 0.034). These results remained significant when controlling for age or sex (Supplementary Text). We additionally examined whether behavioral or brain anatomical measures alone could predict children’s future learning outcomes in response to math interventions and found that our TSI, which combined multimodal information, outperformed unimodal behavioral and brain anatomical measures in predicting intervention-induced math learning (Supplementary Text).   Next, we sought to replicate and demonstrate the generalizability of these findings using a 4-week intervention dataset ( n = 61 children; age: M = 8.2 years, SD = 0.6 years; 33 females; ), which constituted an independent intervention sample that also did not overlap with the cohort used to construct the MAIP . Gray matter volume of 246 brain regions was estimated for each child on the basis of the structural brain image acquired before tutoring. Learning outcome was measured as improvement in the numerical distance effect of processing speed for math problems in a symbolic quantity discrimination task in response to tutoring (Materials and Methods). As described previously, we computed three TSIs in each child and generated predicted scores using a general linear model. As shown in , these three similarity indices for the three components together significantly predicted learning outcomes in the 4-week intervention dataset ( r = 0.38, P = 0.003). Similar results were observed using the similarity index of only the dominant first component, with a large TSI associated with intervention-induced improvement in response time for math problems ( r = −0.25, P = 0.048). These results remained consistent when controlling for age or sex (Supplementary Text). Conventional unimodal prediction analysis yielded similar results as the 8-week intervention cohort (Supplementary Text). These results suggest that our multimodal TSI can reliably predict children’s learning outcomes in both intervention cohorts despite differences in the tutoring protocol, markedly outperforming conventional unimodal approaches.   We investigated the interplay between brain structure, transcriptomic profiles, and individual differences in children’s mathematical abilities and responsiveness to interventions. Our investigation yielded three major findings. First, we uncovered robust associations between gray matter volume and individual differences in children’s mathematical abilities, which were replicated across two independent cohorts. Second, we found that gene expression profiles accounted for 72% of the spatial variance in the association of gray matter volume with children’s mathematical ability. Notably, enrichment analysis revealed that the highest contributors to the gene expression profile overlapped with mathematical ability–related, but not reading- and working memory–related, genes. In addition, we showed that these genes were enriched for neurons, synaptic elements, and potassium channel activity. Third, we developed TSI that measured the extent to which each child’s brain anatomy aligned with brain-wide gene expression profiles linked to mathematical abilities. We found that TSI, measured before intervention, predicted learning outcomes in response to math tutoring. This predictive ability was replicated in two different cohorts of children, despite variations in tutoring content and duration and participant profiles. Together, these findings not only delineate robust brain structural features and transcriptomic signatures that underlie mathematical abilities but also reveal the potential of transcriptome-imaging measures in predicting math learning outcomes. Our findings pave the way for a deeper understanding of the neuroanatomical, transcriptomic factors and molecular processes driving mathematical cognition and learning. Neuroanatomical correlates of individual differences in mathematical ability and identification of a robust MAIP Our first goal was to investigate the associations between gray matter volume and individual differences in children’s mathematical abilities. We implemented a comprehensive multivariate analysis strategy across two large, independent cohorts, marking a substantial advance over previous research that was largely limited to univariate analyses focused on individual brain regions and confined to single cohort studies (refer to table S1 for a summary of previous studies). Multivariate canonical correlation analysis yielded a MAIP, characterized by higher gray matter volume associated with lower mathematical abilities in parietal, ventrotemporal, and PFC regions implicated in mathematical cognition . The validity of the MAIP was confirmed by its significant overlap with functional activation maps drawn from meta-analyses of functional neuroimaging (functional MRI) studies using tasks related to numerical reasoning and mental arithmetic. A close correspondence was observed across various math-related terms and math topic maps, reinforcing the robustness and cognitive relevance of the MAIP. Furthermore, spatial correlations between the MAIP and meta-analytic brain activation maps for math-related terms were significantly higher compared to those between the MAIP and control terms and topic maps, such as reading and writing. These results highlight the functional relevance of the MAIP, revealing that mathematical abilities are underpinned by a distributed network of brain regions. Results challenge the notion that mathematical abilities rely on any single region operating in isolation. This insight represents a substantial advance over previous neuroimaging studies linking mathematical cognition specifically to individual brain regions, such as the PPC. Moreover, our study fills a critical gap by elucidating brain regions that collectively underpin individual differences in children’s mathematical abilities. Our findings align with an emerging systems neuroscience perspective, which emphasizes the role of distributed brain networks in mathematical cognition . Our findings also help to resolve inconsistencies in previous research, which have reported mixed findings, ranging from positive , negative , and even nonsignificant correlations  between gray matter and children’s mathematical abilities. A key distinction of our work is the unprecedented validation of our findings across two independent cohorts. Notably, we used a cross-cohort external validation approach, which is increasingly used and highly recommended for demonstrating the generalizability of brain-behavior findings . In this approach, a prediction model is constructed using data from one cohort, and subsequently, this model is used to make predictions in an entirely independent cohort. In our case, we demonstrated that the brain-behavior model using data from the Stanford discovery cohort could predict the mathematical abilities of children within the NKI-RS validation cohort. Our results provide robust and replicable evidence for an imaging phenotype of mathematical abilities in children, overcoming the limitations of prior studies that lacked replication and predictive analyses across independent cohorts. Our identification of consistent patterns of multivariate brain-behavior association deepens understanding of the neuroanatomical correlates of mathematical abilities. Furthermore, it paves the way for the creation of reliable and objective brain-based biomarkers, offering a powerful tool for the early detection of mathematical disabilities. Transcriptomic factors and molecular processes associated with the MAIP Our next goal was to investigate the transcriptomic factors and related molecular processes associated with mathematical abilities. To accomplish this, we used the Allen human brain-wide transcriptomic atlas, which enabled us to probe the relationship between the MAIP, as characterized above, and gene expression profiles. We found that gene expression profiles accounted for more than 70% of the spatial variability in MAIP. This finding is particularly noteworthy when compared to prior transcriptome-imaging studies, which typically accounted for only 20 to 30% of the variance . Our results stand out by demonstrating a substantially higher percentage of variance explained and suggest that genetic influences contribute significantly to mathematical abilities. Our findings are consistent with 60 to 70% heritability of mathematical abilities reported in twin studies . Four key aspects emerge from the transcriptome-imaging analysis. First, our transcriptomic analysis revealed that low mathematical abilities are associated with high gene expression in distributed brain regions involved in mathematical cognition, including the posterior parietal, prefrontal, and ventrotemporal temporal cortices . This pattern of gene expression coincides with areas displaying increased gray matter volume within the MAIP, suggesting that specific, brain-wide gene expression profiles form the foundation of the MAIP. Critically, in contrast to GWAS, our approach captures the variability in gene expression across diverse brain regions, each contributing to distinct cognitive processes . Through the integration of multimodal brain imaging, gene expression, and cognitive measures, our study provides a more comprehensive and precise understanding of transcriptomic factors and molecular processes that drive mathematical abilities. Second, we leveraged the BrainSpan gene expression atlas  to situate our findings within the broader context of brain development, extending beyond the static snapshot of the adult brain provided by the Allen Human Brain Atlas. The BrainSpan atlas of the developing human brain encompasses extensive gene expression data across various developmental stages. This comprehensive atlas includes the temporal and spatial transcriptional organization of the human brain from 8 weeks after conception to 40 years, allowing us to trace changes in gene expression across development. We plotted the developmental trajectory of gene expression by projecting the gene weights for MAIP derived from partial least squares regression analysis onto the gene expression data from BrainSpan, for each of the six developmental periods. The age range of the late childhood group overlapped with that of our sample. This analysis revealed the stability of gene expression associated with the MAIP from late childhood to adulthood. It suggests that the age differences between the donors in the Allen Human Brain Atlas and our developmental study sample are unlikely to have markedly influenced our findings. This result is broadly consistent with the previous observation of the high heritability of the stability of educational achievement from childhood to adolescence . Third, the genes whose expression profiles most strongly correlated with MAIP were not randomly distributed. Rather, they were significantly enriched by candidate genes previously reported in one of the largest GWAS of cognitive function ( n = 1.1 million individuals), which identified 618 and 365 SNPs associated with mathematical ability and the highest math class taken . Notably, the MAIP genes were not significantly enriched by candidate genes associated with reading or working memory . This result highlights the domain-specific functional relevance of these genes to mathematical cognition. Our findings suggest that both allelic variation and brain-wide expression profiles converge on a specific set of candidate genes linked to mathematical abilities. Future research may further clarify the specific mechanisms by which the genes within this set contribute to mathematical abilities. Fourth, Gene Ontology analysis revealed that the genes whose expression profiles most strongly correlated with MAIP were also significantly enriched for neuron and synapse parts. This enrichment aligns with evidence from studies in mice and nonhuman primates, which have shown that neuronal and synaptic plasticity are critical mechanisms underlying volumetric changes in the brain . These changes often occur in response to experiences and are thought to be foundational in learning and skill acquisition . Furthermore, previous research has shown that synaptic ion channels also modulate communication between neurons, a mechanism essential for learning . In humans, postmortem studies have corroborated this link between brain structure and neuronal and synaptic plasticity . Our findings extend this understanding by identifying a specific pattern of gene expression that aligns with these physiological phenomena, thereby offering further insight into the transcriptomic and molecular underpinnings of the MAIP and, by extension, mathematical abilities. Gene ontology analysis further revealed that the genes whose expression profiles most strongly correlated with MAIP were also enriched in domains associated with potassium channel activity. This aligns with the large GWAS by Lee et al. , which implicated several potassium channel–associated genes in mathematical ability, including KCNN2 , KCNB1 , KCNJ6 , and KCNJ4 . These genes are known to exhibit high expression levels in the brain, influencing the regulation of neuronal excitability . Specific genetic anomalies, such as a frameshift deletion in KCNN2 , have been found in individuals with learning disabilities  and neurodevelopmental delays . Other related studies have pinpointed connections between potassium channelopathies and intellectual disability . These observations, combined with our findings, underscore the critical role of voltage-gated potassium channels in learning, memory, and complex cognitive functions and suggest a transcriptomic and molecular mechanism underlying individual variation in mathematical abilities. Our findings further suggest that altered expression of potassium channel–associated genes might contribute to mathematical disabilities. Considering the high heritability and persistent nature of mathematical disabilities , our study highlights a potential future research avenue in the neuropathology of dyscalculia . Moreover, because noninvasive brain stimulation techniques can alter potassium conductance , our findings may also inspire the development of innovative brain stimulation strategies to remediate learning difficulties, moving beyond traditional behavioral interventions. Together, these findings illuminate the transcriptomic and molecular correlates of the MAIP, enriching our understanding of the biological mechanisms underlying mathematical abilities. Transcriptomic correlates of MAIP predict learning outcomes following math tutoring Our final goal was to determine whether transcriptomic correlates of MAIP affect learning and to assess its prognostic value in forecasting learning outcomes in two independent cohorts of children who participated in math tutoring programs. We developed a TSI to quantify the alignment between MAIP-related gene expression profiles and individual gray matter volume patterns within these cohorts. We then tested the efficacy of TSI in predicting children’s learning outcomes across both cohorts . The TSI metric acts as a crucial link, connecting the MAIP with changes in math performance due to training, where higher TSI values signify greater congruence between math ability–related transcriptomic factors and brain structure. Children in the first intervention participated in an 8-week math tutoring program that aimed to enhance arithmetic problem solving. This involved three weekly 1-hour sessions, encompassing 4 lessons on basic number properties, 2 on counting, and 15 on number families, as well as two review sessions . The second intervention spanned 4 weeks and aimed to enhance number sense. This involved three weekly 1-hour sessions, covering topics such as counting principles, enumeration, and nonsymbolic and symbolic number comparisons . Our analysis revealed that higher TSI consistently predicted greater learning gains in these varied intervention cohorts, despite differences in the content and duration of the tutoring programs, outperforming conventional approaches. Thus, children whose brain structures more closely align with MAIP-related gene expression profiles may have neural substrates that are more conducive to mathematical learning, leading to more effective responses to educational interventions. These results elucidate how gene expression shapes academic performance and learning . Previous genetic studies have revealed high heritability and identified candidate genes associated with mathematical abilities . Our findings go beyond cross-sectional studies and overcome limitations of heritability and GWAS analyses, by providing evidence for the role of gene expression in learning following tutoring. Our study identifies and offers a quantitative framework for reliable biomarkers that predict individual differences in response to interventions, regardless of the type of intervention delivered. Our findings also provide direct evidence for molecular factors that underlie both mathematical abilities and, more importantly, mathematical learning. A potential link between genetic mechanisms of mathematical abilities and learning has been previously observed in an educational context by Harden et al. , which used a polygenic score based on candidate genes previously identified in the GWAS by Lee and colleagues . The researchers found that students with higher polygenic scores in the ninth grade were more likely to pursue advanced mathematical classes by the end of high school. Our study corroborates and extends these findings more directly in the context of short-term math tutoring protocols and identifies the TSI as a robust predictor of responsiveness to interventions designed to remediate mathematical difficulties. The development of reliable transcriptome-based biomarkers may help determine which children are more likely to benefit from these interventions and address the pressing need to identify children who may require more intensive or alternative interventions . More generally, these findings bridge the conceptual and methodological gap between molecular- and macro-level neurobiological mechanisms that underpin mathematical abilities and learning. Our study, while providing critical insights into the neuroanatomical, transcriptomic, and molecular mechanisms of mathematical abilities and learning, has limitations that warrant consideration and point to avenues for future work. We focused on school-age children between 7 and 13 years old, a crucial phase for the acquiring foundational mathematical abilities . Additional research with larger longitudinal datasets is needed to clarify how the neuroanatomical, transcriptomic, and molecular mechanisms of mathematical abilities and learning change with age. Furthermore, Allen Human Brain Atlas donors were older adults. To mitigate this discrepancy, we supplemented our analysis with data from the BrainSpan Developmental Gene Atlas, which included data from children. Although the utilization of BrainSpan revealed relatively stable gene expression from late childhood to adulthood, it is worth noting that this gene expression data has limited coverage of brain regions. Further research with gene expression data spanning the whole brain in the developing brain is needed to systematically investigate how these transcriptomic and molecular mechanisms shape functionally relevant brain structures to support math-related information processing across various stages of development. Currently, such investigation remains a formidable challenge, as constructing such atlases requires access to difficult-to-obtain postmortem brain samples from children and acquiring brain-wide gene expression profiles requires substantial resources and specialized techniques. Our research identifies a significant relation between individual differences in mathematical abilities and the MAIP, independent of the influences of working memory. To more fully dissect how fundamental cognitive operations affect the MAIP, future research will need to use more comprehensive and in-depth cognitive evaluations. These assessments are essential for identifying the precise cognitive mechanisms supporting the MAIP and clarifying the complex interactions between cognitive abilities and mathematical proficiency. We investigated the intricate interplay between neuroanatomical, transcriptomic, and molecular factors in shaping mathematical abilities and learning in school age children. Our multimodal analysis established a replicable and generalizable MAIP and elucidated the molecular underpinnings of mathematical abilities by identifying transcriptomic influences on the MAIP. With gene expression profiles explaining more than 70% of spatial variability in MAIP, our study demonstrated a connection between gene expression and mathematical proficiency. TSI reliably predicted learning outcomes across mathematical intervention contexts and cohorts. Together, our findings advance knowledge of the interplay between neuroanatomical, transcriptomic, and molecular correlates of mathematical ability, bridge micro- and macro-level brain and cognitive functions, and uncover biomarkers of learning, with implications for personalized educational interventions.   Our first goal was to investigate the associations between gray matter volume and individual differences in children’s mathematical abilities. We implemented a comprehensive multivariate analysis strategy across two large, independent cohorts, marking a substantial advance over previous research that was largely limited to univariate analyses focused on individual brain regions and confined to single cohort studies (refer to table S1 for a summary of previous studies). Multivariate canonical correlation analysis yielded a MAIP, characterized by higher gray matter volume associated with lower mathematical abilities in parietal, ventrotemporal, and PFC regions implicated in mathematical cognition . The validity of the MAIP was confirmed by its significant overlap with functional activation maps drawn from meta-analyses of functional neuroimaging (functional MRI) studies using tasks related to numerical reasoning and mental arithmetic. A close correspondence was observed across various math-related terms and math topic maps, reinforcing the robustness and cognitive relevance of the MAIP. Furthermore, spatial correlations between the MAIP and meta-analytic brain activation maps for math-related terms were significantly higher compared to those between the MAIP and control terms and topic maps, such as reading and writing. These results highlight the functional relevance of the MAIP, revealing that mathematical abilities are underpinned by a distributed network of brain regions. Results challenge the notion that mathematical abilities rely on any single region operating in isolation. This insight represents a substantial advance over previous neuroimaging studies linking mathematical cognition specifically to individual brain regions, such as the PPC. Moreover, our study fills a critical gap by elucidating brain regions that collectively underpin individual differences in children’s mathematical abilities. Our findings align with an emerging systems neuroscience perspective, which emphasizes the role of distributed brain networks in mathematical cognition . Our findings also help to resolve inconsistencies in previous research, which have reported mixed findings, ranging from positive , negative , and even nonsignificant correlations  between gray matter and children’s mathematical abilities. A key distinction of our work is the unprecedented validation of our findings across two independent cohorts. Notably, we used a cross-cohort external validation approach, which is increasingly used and highly recommended for demonstrating the generalizability of brain-behavior findings . In this approach, a prediction model is constructed using data from one cohort, and subsequently, this model is used to make predictions in an entirely independent cohort. In our case, we demonstrated that the brain-behavior model using data from the Stanford discovery cohort could predict the mathematical abilities of children within the NKI-RS validation cohort. Our results provide robust and replicable evidence for an imaging phenotype of mathematical abilities in children, overcoming the limitations of prior studies that lacked replication and predictive analyses across independent cohorts. Our identification of consistent patterns of multivariate brain-behavior association deepens understanding of the neuroanatomical correlates of mathematical abilities. Furthermore, it paves the way for the creation of reliable and objective brain-based biomarkers, offering a powerful tool for the early detection of mathematical disabilities.   Our next goal was to investigate the transcriptomic factors and related molecular processes associated with mathematical abilities. To accomplish this, we used the Allen human brain-wide transcriptomic atlas, which enabled us to probe the relationship between the MAIP, as characterized above, and gene expression profiles. We found that gene expression profiles accounted for more than 70% of the spatial variability in MAIP. This finding is particularly noteworthy when compared to prior transcriptome-imaging studies, which typically accounted for only 20 to 30% of the variance . Our results stand out by demonstrating a substantially higher percentage of variance explained and suggest that genetic influences contribute significantly to mathematical abilities. Our findings are consistent with 60 to 70% heritability of mathematical abilities reported in twin studies . Four key aspects emerge from the transcriptome-imaging analysis. First, our transcriptomic analysis revealed that low mathematical abilities are associated with high gene expression in distributed brain regions involved in mathematical cognition, including the posterior parietal, prefrontal, and ventrotemporal temporal cortices . This pattern of gene expression coincides with areas displaying increased gray matter volume within the MAIP, suggesting that specific, brain-wide gene expression profiles form the foundation of the MAIP. Critically, in contrast to GWAS, our approach captures the variability in gene expression across diverse brain regions, each contributing to distinct cognitive processes . Through the integration of multimodal brain imaging, gene expression, and cognitive measures, our study provides a more comprehensive and precise understanding of transcriptomic factors and molecular processes that drive mathematical abilities. Second, we leveraged the BrainSpan gene expression atlas  to situate our findings within the broader context of brain development, extending beyond the static snapshot of the adult brain provided by the Allen Human Brain Atlas. The BrainSpan atlas of the developing human brain encompasses extensive gene expression data across various developmental stages. This comprehensive atlas includes the temporal and spatial transcriptional organization of the human brain from 8 weeks after conception to 40 years, allowing us to trace changes in gene expression across development. We plotted the developmental trajectory of gene expression by projecting the gene weights for MAIP derived from partial least squares regression analysis onto the gene expression data from BrainSpan, for each of the six developmental periods. The age range of the late childhood group overlapped with that of our sample. This analysis revealed the stability of gene expression associated with the MAIP from late childhood to adulthood. It suggests that the age differences between the donors in the Allen Human Brain Atlas and our developmental study sample are unlikely to have markedly influenced our findings. This result is broadly consistent with the previous observation of the high heritability of the stability of educational achievement from childhood to adolescence . Third, the genes whose expression profiles most strongly correlated with MAIP were not randomly distributed. Rather, they were significantly enriched by candidate genes previously reported in one of the largest GWAS of cognitive function ( n = 1.1 million individuals), which identified 618 and 365 SNPs associated with mathematical ability and the highest math class taken . Notably, the MAIP genes were not significantly enriched by candidate genes associated with reading or working memory . This result highlights the domain-specific functional relevance of these genes to mathematical cognition. Our findings suggest that both allelic variation and brain-wide expression profiles converge on a specific set of candidate genes linked to mathematical abilities. Future research may further clarify the specific mechanisms by which the genes within this set contribute to mathematical abilities. Fourth, Gene Ontology analysis revealed that the genes whose expression profiles most strongly correlated with MAIP were also significantly enriched for neuron and synapse parts. This enrichment aligns with evidence from studies in mice and nonhuman primates, which have shown that neuronal and synaptic plasticity are critical mechanisms underlying volumetric changes in the brain . These changes often occur in response to experiences and are thought to be foundational in learning and skill acquisition . Furthermore, previous research has shown that synaptic ion channels also modulate communication between neurons, a mechanism essential for learning . In humans, postmortem studies have corroborated this link between brain structure and neuronal and synaptic plasticity . Our findings extend this understanding by identifying a specific pattern of gene expression that aligns with these physiological phenomena, thereby offering further insight into the transcriptomic and molecular underpinnings of the MAIP and, by extension, mathematical abilities. Gene ontology analysis further revealed that the genes whose expression profiles most strongly correlated with MAIP were also enriched in domains associated with potassium channel activity. This aligns with the large GWAS by Lee et al. , which implicated several potassium channel–associated genes in mathematical ability, including KCNN2 , KCNB1 , KCNJ6 , and KCNJ4 . These genes are known to exhibit high expression levels in the brain, influencing the regulation of neuronal excitability . Specific genetic anomalies, such as a frameshift deletion in KCNN2 , have been found in individuals with learning disabilities  and neurodevelopmental delays . Other related studies have pinpointed connections between potassium channelopathies and intellectual disability . These observations, combined with our findings, underscore the critical role of voltage-gated potassium channels in learning, memory, and complex cognitive functions and suggest a transcriptomic and molecular mechanism underlying individual variation in mathematical abilities. Our findings further suggest that altered expression of potassium channel–associated genes might contribute to mathematical disabilities. Considering the high heritability and persistent nature of mathematical disabilities , our study highlights a potential future research avenue in the neuropathology of dyscalculia . Moreover, because noninvasive brain stimulation techniques can alter potassium conductance , our findings may also inspire the development of innovative brain stimulation strategies to remediate learning difficulties, moving beyond traditional behavioral interventions. Together, these findings illuminate the transcriptomic and molecular correlates of the MAIP, enriching our understanding of the biological mechanisms underlying mathematical abilities.   Our final goal was to determine whether transcriptomic correlates of MAIP affect learning and to assess its prognostic value in forecasting learning outcomes in two independent cohorts of children who participated in math tutoring programs. We developed a TSI to quantify the alignment between MAIP-related gene expression profiles and individual gray matter volume patterns within these cohorts. We then tested the efficacy of TSI in predicting children’s learning outcomes across both cohorts . The TSI metric acts as a crucial link, connecting the MAIP with changes in math performance due to training, where higher TSI values signify greater congruence between math ability–related transcriptomic factors and brain structure. Children in the first intervention participated in an 8-week math tutoring program that aimed to enhance arithmetic problem solving. This involved three weekly 1-hour sessions, encompassing 4 lessons on basic number properties, 2 on counting, and 15 on number families, as well as two review sessions . The second intervention spanned 4 weeks and aimed to enhance number sense. This involved three weekly 1-hour sessions, covering topics such as counting principles, enumeration, and nonsymbolic and symbolic number comparisons . Our analysis revealed that higher TSI consistently predicted greater learning gains in these varied intervention cohorts, despite differences in the content and duration of the tutoring programs, outperforming conventional approaches. Thus, children whose brain structures more closely align with MAIP-related gene expression profiles may have neural substrates that are more conducive to mathematical learning, leading to more effective responses to educational interventions. These results elucidate how gene expression shapes academic performance and learning . Previous genetic studies have revealed high heritability and identified candidate genes associated with mathematical abilities . Our findings go beyond cross-sectional studies and overcome limitations of heritability and GWAS analyses, by providing evidence for the role of gene expression in learning following tutoring. Our study identifies and offers a quantitative framework for reliable biomarkers that predict individual differences in response to interventions, regardless of the type of intervention delivered. Our findings also provide direct evidence for molecular factors that underlie both mathematical abilities and, more importantly, mathematical learning. A potential link between genetic mechanisms of mathematical abilities and learning has been previously observed in an educational context by Harden et al. , which used a polygenic score based on candidate genes previously identified in the GWAS by Lee and colleagues . The researchers found that students with higher polygenic scores in the ninth grade were more likely to pursue advanced mathematical classes by the end of high school. Our study corroborates and extends these findings more directly in the context of short-term math tutoring protocols and identifies the TSI as a robust predictor of responsiveness to interventions designed to remediate mathematical difficulties. The development of reliable transcriptome-based biomarkers may help determine which children are more likely to benefit from these interventions and address the pressing need to identify children who may require more intensive or alternative interventions . More generally, these findings bridge the conceptual and methodological gap between molecular- and macro-level neurobiological mechanisms that underpin mathematical abilities and learning. Our study, while providing critical insights into the neuroanatomical, transcriptomic, and molecular mechanisms of mathematical abilities and learning, has limitations that warrant consideration and point to avenues for future work. We focused on school-age children between 7 and 13 years old, a crucial phase for the acquiring foundational mathematical abilities . Additional research with larger longitudinal datasets is needed to clarify how the neuroanatomical, transcriptomic, and molecular mechanisms of mathematical abilities and learning change with age. Furthermore, Allen Human Brain Atlas donors were older adults. To mitigate this discrepancy, we supplemented our analysis with data from the BrainSpan Developmental Gene Atlas, which included data from children. Although the utilization of BrainSpan revealed relatively stable gene expression from late childhood to adulthood, it is worth noting that this gene expression data has limited coverage of brain regions. Further research with gene expression data spanning the whole brain in the developing brain is needed to systematically investigate how these transcriptomic and molecular mechanisms shape functionally relevant brain structures to support math-related information processing across various stages of development. Currently, such investigation remains a formidable challenge, as constructing such atlases requires access to difficult-to-obtain postmortem brain samples from children and acquiring brain-wide gene expression profiles requires substantial resources and specialized techniques. Our research identifies a significant relation between individual differences in mathematical abilities and the MAIP, independent of the influences of working memory. To more fully dissect how fundamental cognitive operations affect the MAIP, future research will need to use more comprehensive and in-depth cognitive evaluations. These assessments are essential for identifying the precise cognitive mechanisms supporting the MAIP and clarifying the complex interactions between cognitive abilities and mathematical proficiency. We investigated the intricate interplay between neuroanatomical, transcriptomic, and molecular factors in shaping mathematical abilities and learning in school age children. Our multimodal analysis established a replicable and generalizable MAIP and elucidated the molecular underpinnings of mathematical abilities by identifying transcriptomic influences on the MAIP. With gene expression profiles explaining more than 70% of spatial variability in MAIP, our study demonstrated a connection between gene expression and mathematical proficiency. TSI reliably predicted learning outcomes across mathematical intervention contexts and cohorts. Together, our findings advance knowledge of the interplay between neuroanatomical, transcriptomic, and molecular correlates of mathematical ability, bridge micro- and macro-level brain and cognitive functions, and uncover biomarkers of learning, with implications for personalized educational interventions.   Experimental design The study investigated the neurobiological mechanisms that drive individual differences in mathematical abilities and learning, using four datasets of children (ages 7 to 13 years) and a three-part analysis . The first part examined the structural brain organization correlates of individual differences in mathematical abilities, including a discovery analysis in the Stanford cohort and a replication/validation analysis in the NKI-RS cohort. The second part examined gene expression profiles associated with the structural brain organization underlying individual differences in mathematical abilities by using a partial least squares regression analysis, using human brain-wide gene expression data from the Allen Institute for Brain Science. The third part examined whether the coupling between mathematical ability–related gene expression profiles and a child’s structural brain organization can constrain or shape their mathematical learning in two math intervention datasets with distinct tutoring protocols. Participants We used four cohorts of children in this study. Two of these cohorts, namely, the Stanford discovery cohort and NKI-RS validation cohort, were used to identify and replicate the MAIP. The other two cohorts, namely, the 8-week math intervention cohort and 4-week math intervention cohort, were used to investigate the MAIP and the underlying gene expression–based predictor for mathematical learning. Except for the NKI-RS validation cohort, all children were recruited from the San Francisco Bay Area through flyers or poster advertisements across multiple school districts. Informed consent was obtained from the legal guardian of each child. All study protocols were approved by the Stanford University Institutional Review Board. The inclusion criteria included having a normal full-scale intelligence quotient (≥80) as measured by the Wechsler Abbreviated Scale of Intelligence (WASI-I or WASI-II), with no history of claustrophobia, head injury, serious neurological or medical illness, psychosis, mania/bipolar disorder, major depression, substance abuse, sensory impairment, birth weight less than 2000 g, and/or gestational age of less than 34 weeks. The four cohorts are described below, and their demographic and neuropsychological characteristics are summarized in and . Stanford discovery cohort A total of 284 children aged 7 to 13 years completed standardized neuropsychological assessment and structural brain MRI scan sessions. The neuropsychological and brain imaging data were collected in the same visit or two visits within a year, considering the duration considered to be stable standardized scores from the WASI-I or WASI-II and the WIAT-II. Data from 65 children were excluded because they did not meet the inclusion criteria, had administration errors in neuropsychological assessments, had poor structural brain image quality, or had tissue segmentation errors during brain gray matter volume estimation. Ultimately, 219 children were included in the analyses. NKI-RS validation cohort We assembled a validation cohort of children whose age range matched that of the Stanford cohort using the enhanced NKI-RS , a large-scale community sample of participants across the life span. Specifically, we queried for children aged 7 to 13 years who had both structural MRI brain imaging and behavioral data of interest. On the basis of demographic and neuropsychological information available in NKI-RS and following the inclusion criteria described above, 91 children were included in the validation cohort. 8-week math intervention cohort Twenty-four children in grade 3 completed an 8-week tutoring study. All children participated in the neuropsychological assessment session, underwent MRI scanning sessions before and after math tutoring, and engaged in an 8-week math-tutoring program. Tutoring sessions occurred three times per week, each of the tutoring sessions taking approximately 1 hour. The tutoring included 4 lessons on basic number properties, 2 lessons on counting, 15 lessons on number families, and 2 lessons on review. Response to tutoring was examined using an arithmetic verification task that assessed reaction time before and after tutoring. More details about the tutoring material, protocols, and task design are provided in . 4-week math intervention cohort A total of 61 children in grade 2 and grade 3 participated in a 4-week tutoring study. All children participated in the neuropsychological assessment session, underwent MRI scanning sessions before and after math tutoring, and engaged in a 4-week math-tutoring program. This cohort did not overlap with the 8-week intervention cohort. Tutoring sessions occurred three times per week, each of the tutoring sessions taking approximately 1 hour. The tutoring included counting principles, nonsymbolic and symbolic enumeration, and comparisons. Response to tutoring was examined using a symbolic quantity discrimination task acquired before and after training. More details about the tutoring materials, protocols and task design are provided by Chang et al. . Standardized assessments of mathematical abilities To measure children’s mathematical abilities, we used grade-normed standardized scores from two subtests for mathematics in the WIAT-II neuropsychological assessment: (i) numerical operations and (ii) math reasoning. In the numerical operations subtest, each child was required to identify and write numbers as well as solve written calculation problems and simple equations including addition, subtraction, multiplication, and division by writing the answers in the provided response booklet. In the math reasoning subtest, each child was required to identify geometric shapes and solve single- and multistep word problems involving time, money, and measurement, by answering a series of problems presented verbally and visually. Assessments of learning outcomes The math task performed during MRI scanning before and after tutoring was used for learning outcome assessment in two intervention cohorts. We used the changes in processing speed to estimate the learning outcomes, considering potential celling effects on accuracy in the two intervention datasets. For the 8-week intervention cohort, arithmetic verification tasks involving solving two blocks with 24 single-digit addition problems were performed during MRI scanning. For each child, the averaged response time for correct trials was calculated at the pretutoring and posttutoring, respectively. The pretutoring response time was subtracted from the posttutoring response time, yielding an index of response time change related to tutoring. The lower (more negative) value of this index indicated better learning outcomes with tutoring. For the 4-week intervention cohort, a symbolic quantity discrimination task involving 64 trials of comparison of magnitude between two symbolic numbers (1 through 9, excluding 5) was performed during MRI scanning. Half trials had a near distance between two numbers (e.g., 6 and 7), while the other half had a far distance (e.g., 2 and 7). Thus, the numerical distance effect was calculated by subtracting the averaged response time for trials with far distance from those with near distance. For each child, the numerical distance effect at pretutoring was subtracted from the numerical distance effect at posttutoring, yielding an index of change in numerical distance effect in response to tutoring. A lower value for this index indicated better learning outcomes due to tutoring. MRI data acquisition High-resolution T1-weighted spoiled grass gradient recalled inversion recovery three-dimensional MRI sequence was acquired in each child on a 3-T Signa scanner (General Electric) at the Richard M. Lucas Center for Imaging at Stanford University with a custom-built eight-channel head coil for Stanford, 8-week intervention, and 4-week intervention cohorts. Head movement was minimized during the scan by placing cushions around the children’s head. During scanning tasks, children held a custom-made MR-compatible computer mouse or button box in their right hand to respond to math tasks. MPRAGE sequence was acquired in each child on a 3.0 T Siemens TIM Trio at the Nathan Kline Institute with a 32-channel head coil for NKI-RS cohort. More details of NKI imaging protocols can be found at http://fcon_1000.projects.nitrc.org/indi/pro/nki.html . The detailed parameters of four cohorts are summarized in table S6. Brain gray matter volume estimation To estimate the gray matter volume for each child, we used the Computational Anatomy Toolbox ( http://dbm.neuro.uni-jena.de/cat12/ ) for SPM12 ( www.fil.ion.ucl.ac.uk/spm/ ). The voxel-based morphometry analysis was performed with the default parameter setting. Briefly, the structural T1 images were denoised with a spatial-adaptive nonlocal means denoising filter approach and resampled with trilinear interpolation to isotropic resolution of 1 mm by 1 mm by 1 mm, followed by bias correction and affine registration. Then, the resampled images were segmented into gray matter, white matter, and cerebrospinal fluid compartments and spatially normalized to Montreal Neurological Institute (MNI) standard space with the Diffeomorphic Anatomical Registration Through Exponentiated Lie algebra method. The segmented modulated images for gray matter were obtained by multiplying the voxel values with the Jacobian determinant derived from the spatial normalization. In this study, we used the modulated images for gray matter, which represent the absolute gray matter volume rather than gray matter density or the relative scaling volume. Last, the segmented modulated gray matter images were smoothed with an isotropic Gaussian kernel of 12-mm full width at half maximum. To facilitate the association analysis between macroscopic brain structures and microscopic transcriptional features as well as to reduce computational demands, we defined 246 regions of interest by using the Brainnetome Atlas ( https://atlas.brainnetome.org ). For a given region, the regional gray matter volume was obtained by averaging the values of all voxels within that region. Association analysis between brain gray matter volume and mathematical abilities To determine the multivariate association between gray matter volume and mathematical ability across children, we performed a canonical correlation analysis, a widely used method for examining brain-behavior associations . Gray matter volume values from 246 regions of interest and two standardized scores of mathematical abilities (WIAT-II numerical operations and math reasoning scores) for Stanford cohort were included as two input variable sets ( X and Y ). The canonical correlation analysis can determine the modes that relate two input variable sets by identifying the optimal linear combination ( A and B ) in each set that make their correlation ρ maximal ρ = corr ( U , V ) (1) where U and V represent the brain or behavioral canonical variates; the canonical correlation is the correlation ρ between these two canonical variates U = AX (2) V = BY (3) where A and B are canonical vectors to make the linear combination of X and Y to reach the maximal canonical correlation. To avoid overfitting issue of the canonical correlation analysis, we used a principal components analysis on brain variable set before the canonical correlation analysis to perform the dimension reduction. We retained components that explained 80 to 90% of the total variance and use the middle threshold as main result as highly similar results were observed across different threshold (see below). A permutation test was performed to estimate the statistical significance of modes from canonical correlation analysis (1000 times). For each permutation, we randomized the order of WIAT-II numerical operations and math reasoning scores in the behavioral variable set, reperformed the canonical correlation analysis, and recorded the resulting canonical correlation under this randomized brain-behavior correspondence. A null distribution of the canonical correlation was obtained with 1000 permutations. The P value for mode from canonical correlation analysis was obtained by counting the number of the canonical correlation in null distribution that higher than the observed canonical correlation from true brain-behavior correspondence. Behavioral measure weight and brain measure weight of the canonical correlation mode were obtained by correlating behavioral and brain canonical variates respectively against the original behavioral and brain variable sets. To evaluate the robustness of the observed structural brain organization underlying mathematical ability, we examined the effects of different canonical correlation analysis strategies, including permutation inference , covariable regression, and threshold of dimension reduction. We estimated the similarity of the canonical correlation mode in each condition (alternative) with the main finding (original) by calculating their Pearson’s correlation for brain ( U ) and behavioral ( V ) canonical variates, which are the brain and behavioral canonical variates that characterized the canonical correlation mode. The similarity of the canonical correlation mode under different analysis strategies were described in table S2. To examine the robustness within cohorts, we performed a 10-fold cross-validation analysis within Stanford cohort. To further establish the functional specificity of the MAIP pattern, we developed a prediction approach that involved the multiplication of the gray matter volumes of each child with the canonical vector of the observed canonical correlation mode and then correlated it with alternative measures of mathematical proficiency, and two nonmathematical metrics including working memory and socioeconomic status. The detailed results are described in Supplementary Text. Spatial association analysis between MAIP pattern and math-related meta-analytic maps from NeuroSynth To further confirm the cognitive relevance of the observed MAIP pattern from the canonical correlation mode, we estimated its spatial correlation with math-related brain activation meta-maps available from the NeuroSynth meta-analytic database ( www.neurosynth.org ). First, we identified 5 math-related terms in the NeuroSynth database  (of a total of 1334 terms): “symbolic,” “arithmetic,” “subtraction,” “calculation,” and “addition.” For each term, we downloaded the uniformity test map, which identifies brain regions consistently activated in functional MRI studies highly associated with that term. To facilitate the subsequent spatial correlation analysis with the downloaded math-related term maps from NeuroSynth, we remapped the MAIP pattern at the voxel level by assigning the same value for each voxel within the same region. We then calculated the Spearman correlation coefficient across voxels between the brain measure weights within the observed MAIP pattern and the z -scores within the brain activation meta-analytic map of each term to assess their spatial correlation. The brain activation meta-analytic maps were resliced to match the resolution of the gray matter volume image. A permutation test was used to determine whether the correlation coefficient was significantly higher than expected by chance. In each permutation, the nonzero voxels in the meta-analytic map were randomly reassigned, and the correlation coefficient was recalculated. This permutation was repeated 1000 times, yielding a null distribution for each term, while considering the varying size of nonzero voxels in each term. The statistical significance of the actual correlation coefficient was examined on the basis of this null distribution. Since the available terms in the NeuroSynth database did not exclusively derive from studies involving children, we also performed an additional meta-analysis. This meta-analysis exclusively included children studies from the references list in the NeuroSynth database (table S3) to further examine the potential influence on our results using the API Reference analysis tools provided by Neurosynth ( https://neurosynth.readthedocs.io/en/latest/reference.html ). Then, we calculated the Spearman correlation coefficient between the obtained children-only math-related brain activation map and the MAIP pattern. For visualization, we calculated the mean z -score for the brain map of each term in 20 bins of the MAIP to show how the value of the brain map for each term changes with MAIP pattern. The 20 bins of the MAIP were divided on the basis of the gradient brain weights, with each bin including five percentiles of regions. In addition, we downloaded the brain activation meta-analytic maps of two reading-related terms, namely, “reading” and “written,” from the NeuroSynth database for a control analysis. Here, reading ability was chosen as a control to examine the specificity of math ability–related findings as (i) reading and math abilities are academically acquired skills often assessed together in educational settings; (ii) they share domain-general cognitive processes as demonstrated by overlapping brain circuits between these domains ; and (iii) from a learning disability perspective, difficulties in both reading and math are relatively common . To examine the relative cognitive specificity of the MAIP pattern for math-related terms, we calculated the Spearman correlation coefficient across voxels between the MAIP and the brain meta-analytic maps of these two terms and compared these correlations with those of math-related terms using a permutation test (1000 times). Additional validation analysis was performed by using the NeuroSynth “topic” maps associated with “math” and “reading” topics (Supplementary Text). Replication analysis in NKI-RS cohort To examine the replication and generalization of our findings from the Stanford cohort, we developed a prediction approach that uses the observed canonical correlation mode, i.e., the multivariable relation, on brain structures to predict math performance in an independent NKI-RS cohort. Because numerical operations from the WIAT-II was the only available math-related standardized assessment data for the NKI-RS cohort, our prediction approach was only used to predict WIAT-II numerical operations score. First, we estimated the gray matter volume of 246 regions of interest based on the Brainnetome Atlas in NKI-RS cohort for each child using the same pipeline as in the Stanford cohort. We projected the gray matter volume of 246 regions of interest onto the same dimension using principal component coefficient of Stanford cohort X NKI . dm = X NKI PCA . Coeff dm (4) where X NKI is the original gray matter volume matrix of individuals by regions for NKI-RS cohort, and X NKI . dm is the resulting transformed matrix after dimension reduction using principal components analysis as the same as Stanford cohort. PCA.Coeff represents transformation coefficient matrix from principal components analysis, and dm represents the specific threshold of dimension reduction. Here, dm is the middle threshold between 80 and 90% as the main result, which retained components that explained 85% of the total variance. The effect of this threshold on math performance prediction was examined and reported in Supplementary Text and table S2. Then, we generated a brain score related to mathematical ability for each child by multiplying the canonical vector, denoted as A , with the individual gray matter volume pattern U NKI = X NKI . dm A (5) where A is the canonical vector from canonical correlation mode from Stanford cohort. U NKI can be considered as a canonical correlation mode weighted composite brain score that is related to mathematical ability. To generate the predicted WIAT-II numerical operations scores for NKI-RS cohort, we scaled the composite brain scores based on the mean ( m ) and standard deviation (SD) of WIAT-II numerical operations scores in the Stanford cohort Y pred . NKI = U NKI SD + m (6) Last, we calculated the Pearson’s correlation coefficient between the predicted numerical operations scores and the actual numerical operations scores of the NKI cohort. As a control analysis, we also calculated the Pearson’s correlation coefficient between the predicted scores and the actual word reading standardized scores from WIAT-II in the NKI-RS cohort. Williams’s test was used to compare the correlation coefficient between math and reading. In addition to reading ability, we tested additional control measures to further enhance the specificity of math ability–related findings (Supplementary Text). Association analysis between gene expression profiles and MAIP pattern To explore whether the spatial layout of MAIP can be explained by regional gene transcription profiles, we examined the spatial association between the MAIP and the gene expression data from the Allen Human Brain Atlas. Gene expression data and preprocessing Regional microarray expression data were obtained from six postmortem brains (one female, ages 24 to 57, 42 ± 13) provided by the Allen Human Brain Atlas ( https://human.brain-map.org ). Data were processed with the abagen toolbox (version 0.1.3; https://github.com/rmarkello/abagen ) using a 246-region volumetric atlas in MNI space. First, microarray probes were reannotated. Probes not matched to a valid Entrez ID were discarded. Next, probes were filtered on the basis of their expression intensity relative to background noise, such that probes with intensity less than the background in ≥50% of samples across donors were discarded. When multiple probes indexed the expression of the same gene, we selected and used the probe with the most consistent pattern of regional variation across donors (i.e., differential stability). Here, regions correspond to the structural designations provided in the ontology from the Allen Human Brain Atlas. The MNI coordinates of tissue samples were updated to those generated via nonlinear registration using the Advanced Normalization Tools ( https://github.com/chrisfilo/alleninf ). Samples were assigned to brain regions in the provided atlas only if their MNI coordinates were directly within a voxel belonging to a parcel. All tissue samples not assigned to a brain region in the provided atlas were discarded. Interindividual variation was addressed by normalizing tissue sample expression values across genes using a robust sigmoid function. Normalized expression values were rescaled to the unit interval. Gene expression values were then normalized across tissue samples using an identical procedure. Samples assigned to the same brain region were averaged separately for each donor and then across donors, yielding a regional expression matrix (i.e., 246 regions × 15,633 genes). Notably, only the 235 regions covered by gene expression data of Allen Human Brain Atlas and included into the following analysis. Partial least squares regression analysis The partial least squares regression analysis  has been widely used in recent studies  to explore the association between gene expression profiles and the pattern of brain features. This analysis has the advantage of handling situations where the observations (i.e., regions) are fewer than the predictor variables (i.e., genes) and has well-developed scripts and a toolbox that are suitable for the imaging-transcriptomics association analysis . We performed partial least squares regression analysis with the predictor variables being regional gene expression maps and the response variables being the MAIP pattern. The statistical significance of the variance explained by the partial least squares component was determined using a permutation test (1000 times) taking into consideration the spatial autocorrelations in surrogate brain maps . Specifically, these surrogate brain maps were created on the basis of the same value and spatial distance matrix among all brain regions. For each permutation, we performed partial least squares analysis using the same input of gene expression data but with the surrogate brain map of MAIP. We recorded the variance explained by the components in each permutation. The spatial similarity between MAIP and the weighted gene expression profile of each significant partial least squares component was calculated using Pearson’s correlation. For each significant component, a bootstrapping method was used to correct the estimation error of the weight of each gene . Then, we ranked the genes in both descending and ascending sequences according to their corrected weights, which represent the positive or negative contribution to the component, resulting in two ranked gene lists for a given component in two directions. This analysis was performed by adopting the scripts shared by Whitaker et al.  and further customized as needed. Developmental trajectory analysis of gene expression profiles associated with MAIP pattern We performed a developmental trajectory analysis to explore how changes in gene expression related to MAIP with age, particularly interesting in the maturation direction of gene expression during the school-age period. We used BrainSpan, an atlas of the developing human brain, which includes microarray expression data from 42 postmortem brains spanning from 8 weeks after conception to 40 years old (19 female; www.brainspan.org ). Normalized gene expression data of 524 brain tissue samples across 16 cortical and subcortical regions were provided by BrainSpan. In line with our research’s emphasis on the cerebral cortex and postnatal childhood, we excluded brain tissue samples from the cerebellum and early prenatal periods because of their undifferentiated nature. We limited our analysis to genes that overlapped with those used in the partial least squares regression analysis conducted earlier using the Allen Human Brain Atlas. This filtering process resulted in 13,786 genes and 469 samples spanning 15 regions for the subsequent analysis. First, we categorized samples from various ages into six developmental stages of interest: prenatal (<38 weeks), infants (0 to 2 years), early childhood (2 to 6 years), late childhood (6 to 13 years), adolescence (13 to 19 years), and adulthood (>19 years). We selected these age bins on the basis of their broad recognition as key developmental stages, the significant overlap of one age bin (late childhood) with our neuroimaging sample, and the availability of a sufficient number of brain tissue samples within each bin. We ensured that each bin contained at least three subsamples of different ages. Next, we calculated the MAIP-weighted average gene expression for each region and age bin. This involved multiplying the gene weight associated with MAIP by the expression data from BrainSpan and summing the values across all included genes. The gene weight linked to MAIP was derived from the first significant component identified in the previous partial least squares regression analysis. Through these processes, we generated representative expressions associated with MAIP for each region and age bin. After performing these calculations for all regions and age bins, we constructed a developmental trajectory illustrating the profiles of MAIP-related gene expression changes over the course of development. Gene set enrichment analysis on MAIP genes We conducted a gene set enrichment analysis to determine whether the candidate math ability–related genes reported in a previous large-scale GWAS  were overrepresented among the strongly correlated genes identified in our ordered gene lists. Gene set enrichment analysis can detect this effect on either side of the ranked gene list, encompassing genes that are strongly positively or negative correlated. As a result, the descending or ascending sequences of ranked gene lists would yield the same outcome. In this context, we use the descending sequences of the ranked gene lists as the primary results for visualization. Math-related gene set definition The gene sets associated with mathematical ability were selected from the largest large-scale GWAS conducted to date, involving 1.1 million individuals . This study considered two math-related measures, namely, mathematical ability and the highest math class taken. We downloaded the SNP lists of mathematical ability and highest math class taken from their supplementary materials. To map the SNPs to known genes, we performed gene annotation by using SNPnexus, a web-based variant annotation tool ( www.snp-nexus.org/v4/ ) based on the human genome assembly GRCh38/hg38. Only the genes overlapped by the SNPs were used in the subsequent analysis. Consequently, two math-related gene lists were compiled: one for mathematical ability ( n = 386) and another for the highest math class taken ( n = 248). Detailed gene lists can be found in table S5. As part of our control analysis, we used three gene sets related to reading, as identified by recent GWAS  and one gene set related to working memory, a more general cognitive function, as summarized by a review of genome-wise association studies . Gene set enrichment analysis The gene set enrichment analysis was performed using the clusterProfiler package ( https://bioconductor.org/packages/clusterProfiler/ ). For each gene set, an enrichment score representing the level of enrichment was obtained. Then, a normalized enrichment score ( NES ) was generated by comparing the enrichment score with those from permutation tests to correct for the size of the gene set. We also performed gene set enrichment analysis for the gene lists of reading  and working memory  as control analyses. To further examine the role of the math-related gene set on MAIP, we reperformed partial least regression analysis with the predictor variables as regional gene expression maps of these genes only and the response variable as the MAIP pattern. The same permutation test was used for statistical significance . In addition, we conducted the developmental trajectory analysis for the math-related gene set only similarly as described above. Gene Ontology enrichment analysis on MAIP genes We also performed Gene Ontology enrichment analysis on the ranked gene list to identify enriched Gene Ontology terms by using GOrilla ( https://cbl-gorilla.cs.technion.ac.il ). The analysis encompassed biological process, molecular function, and cellular component categories. Significant enrichment was set to Benjamini-Hochberg FDR-corrected q < 0.001. Associations between the TSI and mathematical learning outcomes We used the 8- and 4-week intervention datasets to investigate whether the association between mathematical ability–related gene expression profiles and a child’s structural brain organization could predict their mathematical learning outcomes. A TSI was developed by calculating the correlation between structural brain organization, quantified by regional gray matter volume, and the observed significant mathematical ability–related gene expression profiles. This TSI measures how similar a child’s structural brain organization is to brain-wide gene expression profiles related to mathematical ability. For each child in 8- and 4-week intervention datasets, we first estimated the gray matter volumes based on the structural brain MRI image before tutoring and computed TSI for each significant gene expression profiles. Then, we predicted the learning outcome scores by using a general linear model, with these TSIs as predictor variables while learning outcomes as response variables in two intervention datasets, separately. We estimated the prediction effect for both the first dominant component and three components together. A permutation test was used to determine statistical significance (1000 times). Control analysis for learning outcome predictions To determine the specificity of our findings, we performed several control analyses. First, we calculated a similarity index between a child’s structural brain organization and the MAIP. The MAIP similarity index was calculated as the correlation between the MAIP pattern and gray matter volumes estimated on the basis of the structural brain image acquired before tutoring for each child. We related this MAIP similarity index to learning outcomes using a general linear model. We then correlated this similarity index with their learning outcomes in two intervention datasets. Second, we examined the predictive ability of behavioral performance before the intervention on the learning outcomes. For the 8-week intervention cohort, we used the numerical operations and math reasoning scores from the WIAT-II, collected before the intervention as behavioral measures. For the 4-week intervention cohort, we used the math fluency scores from the Woodcock-Johnson III Tests of Early Cognitive and Academic Development collected before the intervention as behavioral measure. Subsequently, we correlated between these behavioral measures and the learning outcomes to assess their predictive ability for the learning outcomes using a general linear model. Williams’s test was used to compare the predictive ability of the correlation coefficient between these indices and TSI.   The study investigated the neurobiological mechanisms that drive individual differences in mathematical abilities and learning, using four datasets of children (ages 7 to 13 years) and a three-part analysis . The first part examined the structural brain organization correlates of individual differences in mathematical abilities, including a discovery analysis in the Stanford cohort and a replication/validation analysis in the NKI-RS cohort. The second part examined gene expression profiles associated with the structural brain organization underlying individual differences in mathematical abilities by using a partial least squares regression analysis, using human brain-wide gene expression data from the Allen Institute for Brain Science. The third part examined whether the coupling between mathematical ability–related gene expression profiles and a child’s structural brain organization can constrain or shape their mathematical learning in two math intervention datasets with distinct tutoring protocols.   We used four cohorts of children in this study. Two of these cohorts, namely, the Stanford discovery cohort and NKI-RS validation cohort, were used to identify and replicate the MAIP. The other two cohorts, namely, the 8-week math intervention cohort and 4-week math intervention cohort, were used to investigate the MAIP and the underlying gene expression–based predictor for mathematical learning. Except for the NKI-RS validation cohort, all children were recruited from the San Francisco Bay Area through flyers or poster advertisements across multiple school districts. Informed consent was obtained from the legal guardian of each child. All study protocols were approved by the Stanford University Institutional Review Board. The inclusion criteria included having a normal full-scale intelligence quotient (≥80) as measured by the Wechsler Abbreviated Scale of Intelligence (WASI-I or WASI-II), with no history of claustrophobia, head injury, serious neurological or medical illness, psychosis, mania/bipolar disorder, major depression, substance abuse, sensory impairment, birth weight less than 2000 g, and/or gestational age of less than 34 weeks. The four cohorts are described below, and their demographic and neuropsychological characteristics are summarized in and . Stanford discovery cohort A total of 284 children aged 7 to 13 years completed standardized neuropsychological assessment and structural brain MRI scan sessions. The neuropsychological and brain imaging data were collected in the same visit or two visits within a year, considering the duration considered to be stable standardized scores from the WASI-I or WASI-II and the WIAT-II. Data from 65 children were excluded because they did not meet the inclusion criteria, had administration errors in neuropsychological assessments, had poor structural brain image quality, or had tissue segmentation errors during brain gray matter volume estimation. Ultimately, 219 children were included in the analyses. NKI-RS validation cohort We assembled a validation cohort of children whose age range matched that of the Stanford cohort using the enhanced NKI-RS , a large-scale community sample of participants across the life span. Specifically, we queried for children aged 7 to 13 years who had both structural MRI brain imaging and behavioral data of interest. On the basis of demographic and neuropsychological information available in NKI-RS and following the inclusion criteria described above, 91 children were included in the validation cohort. 8-week math intervention cohort Twenty-four children in grade 3 completed an 8-week tutoring study. All children participated in the neuropsychological assessment session, underwent MRI scanning sessions before and after math tutoring, and engaged in an 8-week math-tutoring program. Tutoring sessions occurred three times per week, each of the tutoring sessions taking approximately 1 hour. The tutoring included 4 lessons on basic number properties, 2 lessons on counting, 15 lessons on number families, and 2 lessons on review. Response to tutoring was examined using an arithmetic verification task that assessed reaction time before and after tutoring. More details about the tutoring material, protocols, and task design are provided in . 4-week math intervention cohort A total of 61 children in grade 2 and grade 3 participated in a 4-week tutoring study. All children participated in the neuropsychological assessment session, underwent MRI scanning sessions before and after math tutoring, and engaged in a 4-week math-tutoring program. This cohort did not overlap with the 8-week intervention cohort. Tutoring sessions occurred three times per week, each of the tutoring sessions taking approximately 1 hour. The tutoring included counting principles, nonsymbolic and symbolic enumeration, and comparisons. Response to tutoring was examined using a symbolic quantity discrimination task acquired before and after training. More details about the tutoring materials, protocols and task design are provided by Chang et al. .   A total of 284 children aged 7 to 13 years completed standardized neuropsychological assessment and structural brain MRI scan sessions. The neuropsychological and brain imaging data were collected in the same visit or two visits within a year, considering the duration considered to be stable standardized scores from the WASI-I or WASI-II and the WIAT-II. Data from 65 children were excluded because they did not meet the inclusion criteria, had administration errors in neuropsychological assessments, had poor structural brain image quality, or had tissue segmentation errors during brain gray matter volume estimation. Ultimately, 219 children were included in the analyses.   We assembled a validation cohort of children whose age range matched that of the Stanford cohort using the enhanced NKI-RS , a large-scale community sample of participants across the life span. Specifically, we queried for children aged 7 to 13 years who had both structural MRI brain imaging and behavioral data of interest. On the basis of demographic and neuropsychological information available in NKI-RS and following the inclusion criteria described above, 91 children were included in the validation cohort.   Twenty-four children in grade 3 completed an 8-week tutoring study. All children participated in the neuropsychological assessment session, underwent MRI scanning sessions before and after math tutoring, and engaged in an 8-week math-tutoring program. Tutoring sessions occurred three times per week, each of the tutoring sessions taking approximately 1 hour. The tutoring included 4 lessons on basic number properties, 2 lessons on counting, 15 lessons on number families, and 2 lessons on review. Response to tutoring was examined using an arithmetic verification task that assessed reaction time before and after tutoring. More details about the tutoring material, protocols, and task design are provided in .   A total of 61 children in grade 2 and grade 3 participated in a 4-week tutoring study. All children participated in the neuropsychological assessment session, underwent MRI scanning sessions before and after math tutoring, and engaged in a 4-week math-tutoring program. This cohort did not overlap with the 8-week intervention cohort. Tutoring sessions occurred three times per week, each of the tutoring sessions taking approximately 1 hour. The tutoring included counting principles, nonsymbolic and symbolic enumeration, and comparisons. Response to tutoring was examined using a symbolic quantity discrimination task acquired before and after training. More details about the tutoring materials, protocols and task design are provided by Chang et al. .   To measure children’s mathematical abilities, we used grade-normed standardized scores from two subtests for mathematics in the WIAT-II neuropsychological assessment: (i) numerical operations and (ii) math reasoning. In the numerical operations subtest, each child was required to identify and write numbers as well as solve written calculation problems and simple equations including addition, subtraction, multiplication, and division by writing the answers in the provided response booklet. In the math reasoning subtest, each child was required to identify geometric shapes and solve single- and multistep word problems involving time, money, and measurement, by answering a series of problems presented verbally and visually.   The math task performed during MRI scanning before and after tutoring was used for learning outcome assessment in two intervention cohorts. We used the changes in processing speed to estimate the learning outcomes, considering potential celling effects on accuracy in the two intervention datasets. For the 8-week intervention cohort, arithmetic verification tasks involving solving two blocks with 24 single-digit addition problems were performed during MRI scanning. For each child, the averaged response time for correct trials was calculated at the pretutoring and posttutoring, respectively. The pretutoring response time was subtracted from the posttutoring response time, yielding an index of response time change related to tutoring. The lower (more negative) value of this index indicated better learning outcomes with tutoring. For the 4-week intervention cohort, a symbolic quantity discrimination task involving 64 trials of comparison of magnitude between two symbolic numbers (1 through 9, excluding 5) was performed during MRI scanning. Half trials had a near distance between two numbers (e.g., 6 and 7), while the other half had a far distance (e.g., 2 and 7). Thus, the numerical distance effect was calculated by subtracting the averaged response time for trials with far distance from those with near distance. For each child, the numerical distance effect at pretutoring was subtracted from the numerical distance effect at posttutoring, yielding an index of change in numerical distance effect in response to tutoring. A lower value for this index indicated better learning outcomes due to tutoring.   High-resolution T1-weighted spoiled grass gradient recalled inversion recovery three-dimensional MRI sequence was acquired in each child on a 3-T Signa scanner (General Electric) at the Richard M. Lucas Center for Imaging at Stanford University with a custom-built eight-channel head coil for Stanford, 8-week intervention, and 4-week intervention cohorts. Head movement was minimized during the scan by placing cushions around the children’s head. During scanning tasks, children held a custom-made MR-compatible computer mouse or button box in their right hand to respond to math tasks. MPRAGE sequence was acquired in each child on a 3.0 T Siemens TIM Trio at the Nathan Kline Institute with a 32-channel head coil for NKI-RS cohort. More details of NKI imaging protocols can be found at http://fcon_1000.projects.nitrc.org/indi/pro/nki.html . The detailed parameters of four cohorts are summarized in table S6.   To estimate the gray matter volume for each child, we used the Computational Anatomy Toolbox ( http://dbm.neuro.uni-jena.de/cat12/ ) for SPM12 ( www.fil.ion.ucl.ac.uk/spm/ ). The voxel-based morphometry analysis was performed with the default parameter setting. Briefly, the structural T1 images were denoised with a spatial-adaptive nonlocal means denoising filter approach and resampled with trilinear interpolation to isotropic resolution of 1 mm by 1 mm by 1 mm, followed by bias correction and affine registration. Then, the resampled images were segmented into gray matter, white matter, and cerebrospinal fluid compartments and spatially normalized to Montreal Neurological Institute (MNI) standard space with the Diffeomorphic Anatomical Registration Through Exponentiated Lie algebra method. The segmented modulated images for gray matter were obtained by multiplying the voxel values with the Jacobian determinant derived from the spatial normalization. In this study, we used the modulated images for gray matter, which represent the absolute gray matter volume rather than gray matter density or the relative scaling volume. Last, the segmented modulated gray matter images were smoothed with an isotropic Gaussian kernel of 12-mm full width at half maximum. To facilitate the association analysis between macroscopic brain structures and microscopic transcriptional features as well as to reduce computational demands, we defined 246 regions of interest by using the Brainnetome Atlas ( https://atlas.brainnetome.org ). For a given region, the regional gray matter volume was obtained by averaging the values of all voxels within that region.   To determine the multivariate association between gray matter volume and mathematical ability across children, we performed a canonical correlation analysis, a widely used method for examining brain-behavior associations . Gray matter volume values from 246 regions of interest and two standardized scores of mathematical abilities (WIAT-II numerical operations and math reasoning scores) for Stanford cohort were included as two input variable sets ( X and Y ). The canonical correlation analysis can determine the modes that relate two input variable sets by identifying the optimal linear combination ( A and B ) in each set that make their correlation ρ maximal ρ = corr ( U , V ) (1) where U and V represent the brain or behavioral canonical variates; the canonical correlation is the correlation ρ between these two canonical variates U = AX (2) V = BY (3) where A and B are canonical vectors to make the linear combination of X and Y to reach the maximal canonical correlation. To avoid overfitting issue of the canonical correlation analysis, we used a principal components analysis on brain variable set before the canonical correlation analysis to perform the dimension reduction. We retained components that explained 80 to 90% of the total variance and use the middle threshold as main result as highly similar results were observed across different threshold (see below). A permutation test was performed to estimate the statistical significance of modes from canonical correlation analysis (1000 times). For each permutation, we randomized the order of WIAT-II numerical operations and math reasoning scores in the behavioral variable set, reperformed the canonical correlation analysis, and recorded the resulting canonical correlation under this randomized brain-behavior correspondence. A null distribution of the canonical correlation was obtained with 1000 permutations. The P value for mode from canonical correlation analysis was obtained by counting the number of the canonical correlation in null distribution that higher than the observed canonical correlation from true brain-behavior correspondence. Behavioral measure weight and brain measure weight of the canonical correlation mode were obtained by correlating behavioral and brain canonical variates respectively against the original behavioral and brain variable sets. To evaluate the robustness of the observed structural brain organization underlying mathematical ability, we examined the effects of different canonical correlation analysis strategies, including permutation inference , covariable regression, and threshold of dimension reduction. We estimated the similarity of the canonical correlation mode in each condition (alternative) with the main finding (original) by calculating their Pearson’s correlation for brain ( U ) and behavioral ( V ) canonical variates, which are the brain and behavioral canonical variates that characterized the canonical correlation mode. The similarity of the canonical correlation mode under different analysis strategies were described in table S2. To examine the robustness within cohorts, we performed a 10-fold cross-validation analysis within Stanford cohort. To further establish the functional specificity of the MAIP pattern, we developed a prediction approach that involved the multiplication of the gray matter volumes of each child with the canonical vector of the observed canonical correlation mode and then correlated it with alternative measures of mathematical proficiency, and two nonmathematical metrics including working memory and socioeconomic status. The detailed results are described in Supplementary Text.   To further confirm the cognitive relevance of the observed MAIP pattern from the canonical correlation mode, we estimated its spatial correlation with math-related brain activation meta-maps available from the NeuroSynth meta-analytic database ( www.neurosynth.org ). First, we identified 5 math-related terms in the NeuroSynth database  (of a total of 1334 terms): “symbolic,” “arithmetic,” “subtraction,” “calculation,” and “addition.” For each term, we downloaded the uniformity test map, which identifies brain regions consistently activated in functional MRI studies highly associated with that term. To facilitate the subsequent spatial correlation analysis with the downloaded math-related term maps from NeuroSynth, we remapped the MAIP pattern at the voxel level by assigning the same value for each voxel within the same region. We then calculated the Spearman correlation coefficient across voxels between the brain measure weights within the observed MAIP pattern and the z -scores within the brain activation meta-analytic map of each term to assess their spatial correlation. The brain activation meta-analytic maps were resliced to match the resolution of the gray matter volume image. A permutation test was used to determine whether the correlation coefficient was significantly higher than expected by chance. In each permutation, the nonzero voxels in the meta-analytic map were randomly reassigned, and the correlation coefficient was recalculated. This permutation was repeated 1000 times, yielding a null distribution for each term, while considering the varying size of nonzero voxels in each term. The statistical significance of the actual correlation coefficient was examined on the basis of this null distribution. Since the available terms in the NeuroSynth database did not exclusively derive from studies involving children, we also performed an additional meta-analysis. This meta-analysis exclusively included children studies from the references list in the NeuroSynth database (table S3) to further examine the potential influence on our results using the API Reference analysis tools provided by Neurosynth ( https://neurosynth.readthedocs.io/en/latest/reference.html ). Then, we calculated the Spearman correlation coefficient between the obtained children-only math-related brain activation map and the MAIP pattern. For visualization, we calculated the mean z -score for the brain map of each term in 20 bins of the MAIP to show how the value of the brain map for each term changes with MAIP pattern. The 20 bins of the MAIP were divided on the basis of the gradient brain weights, with each bin including five percentiles of regions. In addition, we downloaded the brain activation meta-analytic maps of two reading-related terms, namely, “reading” and “written,” from the NeuroSynth database for a control analysis. Here, reading ability was chosen as a control to examine the specificity of math ability–related findings as (i) reading and math abilities are academically acquired skills often assessed together in educational settings; (ii) they share domain-general cognitive processes as demonstrated by overlapping brain circuits between these domains ; and (iii) from a learning disability perspective, difficulties in both reading and math are relatively common . To examine the relative cognitive specificity of the MAIP pattern for math-related terms, we calculated the Spearman correlation coefficient across voxels between the MAIP and the brain meta-analytic maps of these two terms and compared these correlations with those of math-related terms using a permutation test (1000 times). Additional validation analysis was performed by using the NeuroSynth “topic” maps associated with “math” and “reading” topics (Supplementary Text).   To examine the replication and generalization of our findings from the Stanford cohort, we developed a prediction approach that uses the observed canonical correlation mode, i.e., the multivariable relation, on brain structures to predict math performance in an independent NKI-RS cohort. Because numerical operations from the WIAT-II was the only available math-related standardized assessment data for the NKI-RS cohort, our prediction approach was only used to predict WIAT-II numerical operations score. First, we estimated the gray matter volume of 246 regions of interest based on the Brainnetome Atlas in NKI-RS cohort for each child using the same pipeline as in the Stanford cohort. We projected the gray matter volume of 246 regions of interest onto the same dimension using principal component coefficient of Stanford cohort X NKI . dm = X NKI PCA . Coeff dm (4) where X NKI is the original gray matter volume matrix of individuals by regions for NKI-RS cohort, and X NKI . dm is the resulting transformed matrix after dimension reduction using principal components analysis as the same as Stanford cohort. PCA.Coeff represents transformation coefficient matrix from principal components analysis, and dm represents the specific threshold of dimension reduction. Here, dm is the middle threshold between 80 and 90% as the main result, which retained components that explained 85% of the total variance. The effect of this threshold on math performance prediction was examined and reported in Supplementary Text and table S2. Then, we generated a brain score related to mathematical ability for each child by multiplying the canonical vector, denoted as A , with the individual gray matter volume pattern U NKI = X NKI . dm A (5) where A is the canonical vector from canonical correlation mode from Stanford cohort. U NKI can be considered as a canonical correlation mode weighted composite brain score that is related to mathematical ability. To generate the predicted WIAT-II numerical operations scores for NKI-RS cohort, we scaled the composite brain scores based on the mean ( m ) and standard deviation (SD) of WIAT-II numerical operations scores in the Stanford cohort Y pred . NKI = U NKI SD + m (6) Last, we calculated the Pearson’s correlation coefficient between the predicted numerical operations scores and the actual numerical operations scores of the NKI cohort. As a control analysis, we also calculated the Pearson’s correlation coefficient between the predicted scores and the actual word reading standardized scores from WIAT-II in the NKI-RS cohort. Williams’s test was used to compare the correlation coefficient between math and reading. In addition to reading ability, we tested additional control measures to further enhance the specificity of math ability–related findings (Supplementary Text).   To explore whether the spatial layout of MAIP can be explained by regional gene transcription profiles, we examined the spatial association between the MAIP and the gene expression data from the Allen Human Brain Atlas. Gene expression data and preprocessing Regional microarray expression data were obtained from six postmortem brains (one female, ages 24 to 57, 42 ± 13) provided by the Allen Human Brain Atlas ( https://human.brain-map.org ). Data were processed with the abagen toolbox (version 0.1.3; https://github.com/rmarkello/abagen ) using a 246-region volumetric atlas in MNI space. First, microarray probes were reannotated. Probes not matched to a valid Entrez ID were discarded. Next, probes were filtered on the basis of their expression intensity relative to background noise, such that probes with intensity less than the background in ≥50% of samples across donors were discarded. When multiple probes indexed the expression of the same gene, we selected and used the probe with the most consistent pattern of regional variation across donors (i.e., differential stability). Here, regions correspond to the structural designations provided in the ontology from the Allen Human Brain Atlas. The MNI coordinates of tissue samples were updated to those generated via nonlinear registration using the Advanced Normalization Tools ( https://github.com/chrisfilo/alleninf ). Samples were assigned to brain regions in the provided atlas only if their MNI coordinates were directly within a voxel belonging to a parcel. All tissue samples not assigned to a brain region in the provided atlas were discarded. Interindividual variation was addressed by normalizing tissue sample expression values across genes using a robust sigmoid function. Normalized expression values were rescaled to the unit interval. Gene expression values were then normalized across tissue samples using an identical procedure. Samples assigned to the same brain region were averaged separately for each donor and then across donors, yielding a regional expression matrix (i.e., 246 regions × 15,633 genes). Notably, only the 235 regions covered by gene expression data of Allen Human Brain Atlas and included into the following analysis. Partial least squares regression analysis The partial least squares regression analysis  has been widely used in recent studies  to explore the association between gene expression profiles and the pattern of brain features. This analysis has the advantage of handling situations where the observations (i.e., regions) are fewer than the predictor variables (i.e., genes) and has well-developed scripts and a toolbox that are suitable for the imaging-transcriptomics association analysis . We performed partial least squares regression analysis with the predictor variables being regional gene expression maps and the response variables being the MAIP pattern. The statistical significance of the variance explained by the partial least squares component was determined using a permutation test (1000 times) taking into consideration the spatial autocorrelations in surrogate brain maps . Specifically, these surrogate brain maps were created on the basis of the same value and spatial distance matrix among all brain regions. For each permutation, we performed partial least squares analysis using the same input of gene expression data but with the surrogate brain map of MAIP. We recorded the variance explained by the components in each permutation. The spatial similarity between MAIP and the weighted gene expression profile of each significant partial least squares component was calculated using Pearson’s correlation. For each significant component, a bootstrapping method was used to correct the estimation error of the weight of each gene . Then, we ranked the genes in both descending and ascending sequences according to their corrected weights, which represent the positive or negative contribution to the component, resulting in two ranked gene lists for a given component in two directions. This analysis was performed by adopting the scripts shared by Whitaker et al.  and further customized as needed.   Regional microarray expression data were obtained from six postmortem brains (one female, ages 24 to 57, 42 ± 13) provided by the Allen Human Brain Atlas ( https://human.brain-map.org ). Data were processed with the abagen toolbox (version 0.1.3; https://github.com/rmarkello/abagen ) using a 246-region volumetric atlas in MNI space. First, microarray probes were reannotated. Probes not matched to a valid Entrez ID were discarded. Next, probes were filtered on the basis of their expression intensity relative to background noise, such that probes with intensity less than the background in ≥50% of samples across donors were discarded. When multiple probes indexed the expression of the same gene, we selected and used the probe with the most consistent pattern of regional variation across donors (i.e., differential stability). Here, regions correspond to the structural designations provided in the ontology from the Allen Human Brain Atlas. The MNI coordinates of tissue samples were updated to those generated via nonlinear registration using the Advanced Normalization Tools ( https://github.com/chrisfilo/alleninf ). Samples were assigned to brain regions in the provided atlas only if their MNI coordinates were directly within a voxel belonging to a parcel. All tissue samples not assigned to a brain region in the provided atlas were discarded. Interindividual variation was addressed by normalizing tissue sample expression values across genes using a robust sigmoid function. Normalized expression values were rescaled to the unit interval. Gene expression values were then normalized across tissue samples using an identical procedure. Samples assigned to the same brain region were averaged separately for each donor and then across donors, yielding a regional expression matrix (i.e., 246 regions × 15,633 genes). Notably, only the 235 regions covered by gene expression data of Allen Human Brain Atlas and included into the following analysis.   The partial least squares regression analysis  has been widely used in recent studies  to explore the association between gene expression profiles and the pattern of brain features. This analysis has the advantage of handling situations where the observations (i.e., regions) are fewer than the predictor variables (i.e., genes) and has well-developed scripts and a toolbox that are suitable for the imaging-transcriptomics association analysis . We performed partial least squares regression analysis with the predictor variables being regional gene expression maps and the response variables being the MAIP pattern. The statistical significance of the variance explained by the partial least squares component was determined using a permutation test (1000 times) taking into consideration the spatial autocorrelations in surrogate brain maps . Specifically, these surrogate brain maps were created on the basis of the same value and spatial distance matrix among all brain regions. For each permutation, we performed partial least squares analysis using the same input of gene expression data but with the surrogate brain map of MAIP. We recorded the variance explained by the components in each permutation. The spatial similarity between MAIP and the weighted gene expression profile of each significant partial least squares component was calculated using Pearson’s correlation. For each significant component, a bootstrapping method was used to correct the estimation error of the weight of each gene . Then, we ranked the genes in both descending and ascending sequences according to their corrected weights, which represent the positive or negative contribution to the component, resulting in two ranked gene lists for a given component in two directions. This analysis was performed by adopting the scripts shared by Whitaker et al.  and further customized as needed.   We performed a developmental trajectory analysis to explore how changes in gene expression related to MAIP with age, particularly interesting in the maturation direction of gene expression during the school-age period. We used BrainSpan, an atlas of the developing human brain, which includes microarray expression data from 42 postmortem brains spanning from 8 weeks after conception to 40 years old (19 female; www.brainspan.org ). Normalized gene expression data of 524 brain tissue samples across 16 cortical and subcortical regions were provided by BrainSpan. In line with our research’s emphasis on the cerebral cortex and postnatal childhood, we excluded brain tissue samples from the cerebellum and early prenatal periods because of their undifferentiated nature. We limited our analysis to genes that overlapped with those used in the partial least squares regression analysis conducted earlier using the Allen Human Brain Atlas. This filtering process resulted in 13,786 genes and 469 samples spanning 15 regions for the subsequent analysis. First, we categorized samples from various ages into six developmental stages of interest: prenatal (<38 weeks), infants (0 to 2 years), early childhood (2 to 6 years), late childhood (6 to 13 years), adolescence (13 to 19 years), and adulthood (>19 years). We selected these age bins on the basis of their broad recognition as key developmental stages, the significant overlap of one age bin (late childhood) with our neuroimaging sample, and the availability of a sufficient number of brain tissue samples within each bin. We ensured that each bin contained at least three subsamples of different ages. Next, we calculated the MAIP-weighted average gene expression for each region and age bin. This involved multiplying the gene weight associated with MAIP by the expression data from BrainSpan and summing the values across all included genes. The gene weight linked to MAIP was derived from the first significant component identified in the previous partial least squares regression analysis. Through these processes, we generated representative expressions associated with MAIP for each region and age bin. After performing these calculations for all regions and age bins, we constructed a developmental trajectory illustrating the profiles of MAIP-related gene expression changes over the course of development.   We conducted a gene set enrichment analysis to determine whether the candidate math ability–related genes reported in a previous large-scale GWAS  were overrepresented among the strongly correlated genes identified in our ordered gene lists. Gene set enrichment analysis can detect this effect on either side of the ranked gene list, encompassing genes that are strongly positively or negative correlated. As a result, the descending or ascending sequences of ranked gene lists would yield the same outcome. In this context, we use the descending sequences of the ranked gene lists as the primary results for visualization. Math-related gene set definition The gene sets associated with mathematical ability were selected from the largest large-scale GWAS conducted to date, involving 1.1 million individuals . This study considered two math-related measures, namely, mathematical ability and the highest math class taken. We downloaded the SNP lists of mathematical ability and highest math class taken from their supplementary materials. To map the SNPs to known genes, we performed gene annotation by using SNPnexus, a web-based variant annotation tool ( www.snp-nexus.org/v4/ ) based on the human genome assembly GRCh38/hg38. Only the genes overlapped by the SNPs were used in the subsequent analysis. Consequently, two math-related gene lists were compiled: one for mathematical ability ( n = 386) and another for the highest math class taken ( n = 248). Detailed gene lists can be found in table S5. As part of our control analysis, we used three gene sets related to reading, as identified by recent GWAS  and one gene set related to working memory, a more general cognitive function, as summarized by a review of genome-wise association studies . Gene set enrichment analysis The gene set enrichment analysis was performed using the clusterProfiler package ( https://bioconductor.org/packages/clusterProfiler/ ). For each gene set, an enrichment score representing the level of enrichment was obtained. Then, a normalized enrichment score ( NES ) was generated by comparing the enrichment score with those from permutation tests to correct for the size of the gene set. We also performed gene set enrichment analysis for the gene lists of reading  and working memory  as control analyses. To further examine the role of the math-related gene set on MAIP, we reperformed partial least regression analysis with the predictor variables as regional gene expression maps of these genes only and the response variable as the MAIP pattern. The same permutation test was used for statistical significance . In addition, we conducted the developmental trajectory analysis for the math-related gene set only similarly as described above.   The gene sets associated with mathematical ability were selected from the largest large-scale GWAS conducted to date, involving 1.1 million individuals . This study considered two math-related measures, namely, mathematical ability and the highest math class taken. We downloaded the SNP lists of mathematical ability and highest math class taken from their supplementary materials. To map the SNPs to known genes, we performed gene annotation by using SNPnexus, a web-based variant annotation tool ( www.snp-nexus.org/v4/ ) based on the human genome assembly GRCh38/hg38. Only the genes overlapped by the SNPs were used in the subsequent analysis. Consequently, two math-related gene lists were compiled: one for mathematical ability ( n = 386) and another for the highest math class taken ( n = 248). Detailed gene lists can be found in table S5. As part of our control analysis, we used three gene sets related to reading, as identified by recent GWAS  and one gene set related to working memory, a more general cognitive function, as summarized by a review of genome-wise association studies .   The gene set enrichment analysis was performed using the clusterProfiler package ( https://bioconductor.org/packages/clusterProfiler/ ). For each gene set, an enrichment score representing the level of enrichment was obtained. Then, a normalized enrichment score ( NES ) was generated by comparing the enrichment score with those from permutation tests to correct for the size of the gene set. We also performed gene set enrichment analysis for the gene lists of reading  and working memory  as control analyses. To further examine the role of the math-related gene set on MAIP, we reperformed partial least regression analysis with the predictor variables as regional gene expression maps of these genes only and the response variable as the MAIP pattern. The same permutation test was used for statistical significance . In addition, we conducted the developmental trajectory analysis for the math-related gene set only similarly as described above.   We also performed Gene Ontology enrichment analysis on the ranked gene list to identify enriched Gene Ontology terms by using GOrilla ( https://cbl-gorilla.cs.technion.ac.il ). The analysis encompassed biological process, molecular function, and cellular component categories. Significant enrichment was set to Benjamini-Hochberg FDR-corrected q < 0.001.   We used the 8- and 4-week intervention datasets to investigate whether the association between mathematical ability–related gene expression profiles and a child’s structural brain organization could predict their mathematical learning outcomes. A TSI was developed by calculating the correlation between structural brain organization, quantified by regional gray matter volume, and the observed significant mathematical ability–related gene expression profiles. This TSI measures how similar a child’s structural brain organization is to brain-wide gene expression profiles related to mathematical ability. For each child in 8- and 4-week intervention datasets, we first estimated the gray matter volumes based on the structural brain MRI image before tutoring and computed TSI for each significant gene expression profiles. Then, we predicted the learning outcome scores by using a general linear model, with these TSIs as predictor variables while learning outcomes as response variables in two intervention datasets, separately. We estimated the prediction effect for both the first dominant component and three components together. A permutation test was used to determine statistical significance (1000 times).   To determine the specificity of our findings, we performed several control analyses. First, we calculated a similarity index between a child’s structural brain organization and the MAIP. The MAIP similarity index was calculated as the correlation between the MAIP pattern and gray matter volumes estimated on the basis of the structural brain image acquired before tutoring for each child. We related this MAIP similarity index to learning outcomes using a general linear model. We then correlated this similarity index with their learning outcomes in two intervention datasets. Second, we examined the predictive ability of behavioral performance before the intervention on the learning outcomes. For the 8-week intervention cohort, we used the numerical operations and math reasoning scores from the WIAT-II, collected before the intervention as behavioral measures. For the 4-week intervention cohort, we used the math fluency scores from the Woodcock-Johnson III Tests of Early Cognitive and Academic Development collected before the intervention as behavioral measure. Subsequently, we correlated between these behavioral measures and the learning outcomes to assess their predictive ability for the learning outcomes using a general linear model. Williams’s test was used to compare the predictive ability of the correlation coefficient between these indices and TSI.   