Ph.D Junhao Wen
email: junhaowe@usc.edu

Ilya M Nasrallah

Ahmed Abdulkadir

Theodore D Satterthwaite

Zhijian Yang

Guray Erus

Timothy Robert-Fitzgerald

Ashish Singh

Aristeidis Sotiras

Aleix Boquet-Pujadas

Elizabeth Mamourian

Jimit Doshi

Yuhan Cui

Dhivya Srinivasan

Ioanna Skampardoni

Jiong Chen

Gyujoon Hwang

Mark Bergman

Jingxuan Bao

Yogasudha Veturi

Zhen Zhou

Shu Yang

Paola Dazzan

Rene S Kahn

Hugo G Schnack

Marcus V Zanetti

Eva Meisenzahl

Geraldo F Busatto

Benedicto Crespo-Facorro

Christos Pantelis

Stephen J Wood

Chuanjun Zhuo

Russell T Shinohara

Ruben C Gur

Raquel E Gur

Nikolaos Koutsouleris

Daniel H Wolf

Andrew J Saykin

Marylyn D Ritchie

Li Shen

Paul M Thompson

Olivier Colliot

Katharina Wittfeld

Hans J Grabe

Duygu Tosun

Murat Bilgel

Yang An

Daniel S Marcus

Pamela Lamontagne

Susan R Heckbert

Thomas R Austin

Lenore J Launer

Mark Espeland

Colin L Masters

Paul Maruff

Jurgen Fripp

Sterling C Johnson

John C Morris

Marilyn S Albert

R Nick Bryan

Susan M Resnick

Yong Fan

Mohamad Habes

David Wolk

Haochang Shou

Ph.D Christos Davatzikos

Genomic loci influence patterns of structural covariance in the human brain

Keywords: Pattern of structural covariance, brain imaging genetics, matrix factorization

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Brain structure and function are interrelated via complex networks that operate at multiple scales, ranging from cellular and synaptic processes, such as neural migration, synapse formation, and axon development, to local and broadly connected circuits. [START_REF] Alexander-Bloch | Imaging structural co-variance between human brain regions[END_REF] Due to a fundamental relationship between activity and structure, many normal and pathologic neurobiological processes, driven by genetic and environmental factors, collectively cause coordinated changes in brain morphology.

Structural covariance analyses investigate such coordinated changes by seeking patterns of structural covariation (PSC) across brain regions and individuals. [START_REF] Alexander-Bloch | Imaging structural co-variance between human brain regions[END_REF] For example, during adolescence, PSCs derived from magnetic resonance imaging (MRI) have been considered to reflect a coordinated cortical remodeling as the brain establishes mature networks of functional specialization. [START_REF] Sotiras | Patterns of coordinated cortical remodeling during adolescence and their associations with functional specialization and evolutionary expansion[END_REF] Structural covariance is not only related to normal brain development or aging processes but can also reflect coordinated brain change due to disease. For example, individuals with motor speech dysfunction may develop brain atrophy in Broca's inferior frontal cortex and co-occurring brain atrophy in Wernicke's area of the superior temporal cortex. [START_REF] Blank | Speech production: Wernicke, Broca and beyond[END_REF] Refer to Fig. 1C for an illustrative depiction.

The human brain develops, matures, and degenerates in coordinated patterns of structural covariance at the macrostructural level of brain morphology. [START_REF] Alexander-Bloch | Imaging structural co-variance between human brain regions[END_REF] However, the mechanisms underlying structural covariance are still unclear, and their genetic underpinnings are largely unknown. We hypothesized that brain morphology was driven by multiple genes (i.e., polygenic) collectively operating on different brain areas (i.e., pleiotropic), resulting in connected networks covaried by normal aging and various disease-related processes. Along the causal pathway from underlying genetics to brain morphological changes, we sought to elucidate which genetic underpinnings (e.g., genes), biological processes (e.g., neurogenesis), cellular components (e.g., nuclear membrane), molecular functions (e.g., nucleic acid binding), and neuropathological processes (e.g., Alzheimer's disease) might influence the formation, development, and changes of structural covariance patterns in the human brain.

Previous neuroimaging genome-wide association studies (GWAS) [START_REF] Zhao | Genome-wide association analysis of 19,629 individuals identifies variants influencing regional brain volumes and refines their genetic co-architecture with cognitive and mental health traits[END_REF][START_REF] Elliott | Genome-wide association studies of brain imaging phenotypes in UK Biobank[END_REF] have partially investigated the abovementioned questions and expanded our understanding of the genetic architecture of the human brain. However, they focused on conventional neuroanatomical regions of interest (ROI) instead of data-driven PSCs. In brain imaging research, prior studies have applied structural covariance analysis to elucidate underlying coordinated morphological changes in brain aging and various brain diseases, 1 but have had several limitations. They often relied on pre-defined neuroanatomical ROIs to construct inter-and intra-individual structural covariance networks. These a priori ROIs might not optimally reflect the molecular-functional characteristics of the brain. In addition, most population-based studies have investigated brain structural covariance within a relatively limited scope, such as within relatively small samples, over a relatively narrow age window (e.g., adolescence 2 ), within a single disease (e.g., Parkinson's disease [START_REF] Vignando | Mapping brain structural differences and neuroreceptor correlates in Parkinson's disease visual hallucinations[END_REF] ), or within datasets lacking sufficient diversity in cohort characteristics or MRI scanner protocols. These have been imposed, in part, by limitations in both available cohort size and in the algorithmic implementation of structural covariance analysis, which has been computationally restricted to modest sample sizes when investigated at full image resolution.

Lastly, prior studies have examined brain structural covariance at a single fixed ROI resolution/scale/granularity. While the optimal scale is unknown and may differ by the question of interest, the highly complex organization of the human brain may demonstrate structural covariance patterns that span multiple scales. [START_REF] Bassett | Multiscale Network Organization in the Human Brain[END_REF][START_REF] Schaefer | Local-Global Parcellation of the Human Cerebral Cortex from Intrinsic Functional Connectivity MRI[END_REF] To address this gap, we modified our previously proposed orthogonally projective nonnegative matrix factorization (opNMF [START_REF] Sotiras | Finding imaging patterns of structural covariance via Non-Negative Matrix Factorization[END_REF] ) to its stochastic counterpart, sopNMF. This adaptation allowed us to train the model iteratively on large-scale neuroimaging datasets with a pre-defined number of PSCs (C). Non-negative matrix factorization has gained significant attention in neuroimaging due to its ability to reduce complex data into a sparse, part-based brain representation by projection onto a relatively small number of components (the PSCs). NMF has been shown to substantially improve interpretability and reproducibility compared to other unsupervised methods, such as PCA and ICA, thanks to the non-negative constraint that produces parcellation-like decompositions of complex signals. Our opNMF/sopNMF approach imposed an additional orthonormality constraint 9 (Equation 1 in Method 1), further enhancing sparsity and facilitating clinical interpretability. In our previous work, we applied the opNMF method to 934 youths ages 8-20 to depict the coordinated growth of structural brain networks during adolescence -a period characterized by extensive remodeling of the human cortex to accommodate the rapid expansion of the behavioral repertoire 2 . Remarkably, this study revealed PSCs that exhibited a cortical organization closely aligned with established functional brain networks, such as the well-known 7-network functional parcellation proposed by Yeo et al [START_REF] Yeo | The organization of the human cerebral cortex estimated by intrinsic functional connectivity[END_REF] .

Notably, this alignment emerged without prior assumptions, was data-driven and hypothesisfree, and potentially reflected underlying neurobiological processes related to brain development and aging. Herein, we used large-scale neuroimaging data to investigate the underlying genetic determinant influencing such changes in structural covariance patterns in the human brain.

We examined structural covariance of regional cortical and subcortical volume in the human brain using MRI from a diverse population of 50,699 people from 12 studies, 130 sites, and 12 countries, comprised of cognitively healthy individuals, as well as participants with various diseases/conditions over their lifespan (ages 5 through 97). Herein we present results from coarse to fine scales corresponding to C = 32, 64, 128, 256, 512, and 1024. We hypothesized that PSCs at multiple scales could delineate the human brain's multi-factorial and multi-faceted morphological landscape and genetic architecture in healthy and diseased individuals. We examined the associations between these multi-scale PSCs and common genetic variants at different levels (N=8,469,833 SNPs). In total, 617 novel genomic loci were identified; key pathways (e.g., neurogenesis and reelin signaling) contributed to shaping structural covariance patterns in the human brain. In addition, we leveraged PSCs at multiple scales to better derive individualized imaging signatures of several diseases than any single-scale PSCs using support vector machines. All experimental results and the multi-scale PSCs were integrated into the MuSIC (Multi-scale Structural Imaging Covariance) atlas and made publicly accessible through the BRIDGEPORT (BRaIn knowleDGE PORTal) web portal: https://www.cbica.upenn.edu/bridgeport/.

Results

We summarize this work in three units (I to III) outlined in Fig. 1. In Unit I (Fig. 1A), we present the stochastic orthogonally projective non-negative matrix factorization (sopNMF) algorithm (Method 1), optimized for large-scale multivariate structural covariance analysis. The sopNMF algorithm decomposes large-scale imaging data through online learning to overcome the memory limitations of opNMF. A subgroup of participants with multiple disease diagnoses and healthy controls (ages 5-97, training population, N=4000, Method 2) were sampled from the discovery set (N=32,440, Method 2); their MRI underwent a standard imaging processing pipeline (Method 3A). The processed images were then fit to sopNMF to derive the multi-scale PSCs (N=2003) from the loadings of the factorization (Method 1). We incorporate participants with various disease conditions because previous studies have demonstrated that inter-regional correlated patterns (i.e., depicting a network) show variations in healthy and diseased populations, albeit to a differing degree. that sought to discover associations of PSCs at single nucleotide polymorphism (SNP), gene, or gene set-level; and ii) pattern analysis via support vector machine (Method 5) to derive individualized imaging signatures of several brain diseases and conditions. Unit III (Fig. 1C)

presents BRIDGEPORT, making these massive analytic resources publicly available to the imaging, genomics, and machine learning communities.

Patterns of structural covariance via stochastic orthogonally projective non-negative matrix factorization

We first validated the sopNMF algorithm by showing that it converged to the global minimum of the factorization problem using the comparison population (N=800, Method 2). The sopNMF algorithm achieved similar reconstruction loss and sparsity as opNMF but at reduced memory demand (SI eFigure 1). We then extracted the multi-scale PSCs in the discovery set (N=32,440) and the replication set (N=18,259, Method 2) for Unit II. These PSCs succinctly capture underlying neurobiological processes across the lifespan, including the effects of typical aging processes and various brain diseases. In addition, the multi-scale representation constructs a hierarchy of brain structure networks (e.g., PSCs in cerebellum regions), which models the human brain in a multi-scale topology. [START_REF] Bassett | Multiscale Network Organization in the Human Brain[END_REF][START_REF] Betzel | Multi-scale brain networks[END_REF] Patterns of structural covariance are highly heritable

The multi-scale PSCs are highly heritable (0.05< h 2 <0.78), showing high SNP-based heritability estimates (h 2 ) (Method 4B) for the discovery set (Fig. 2). Specifically, the h 2 estimate was 0.49±0.10, 0.39±0. are presented in BRIDGEPORT. In addition to the abovementioned replication analyses, we also performed several sensitivity analyses (SI eFigure 4a). Our findings revealed the robustness of GWAS signals across both the discovery and replication sets, even when considering four additional brain-related covariates. However, the generalizability of these signals was limited in non-European ancestry populations and independent disease-specific populations (SI eText 1

and SI eFigure 4).

Gene set enrichment analysis highlights pathways that shape patterns of structural covariance

For gene-level associations (Method 4D), we discovered that 164 genes had 2489 gene-PSC pairwise associations with 445 PSCs after Bonferroni correction for the number of genes and PSCs (p-value threshold: 8.6 > -log10[p-value] > 7.1) (SI eFile 9).

Based on these gene-level p-values, we performed hypothesis-free gene set pathway analysis using MAGMA 17 (Method 4E): a more stringent correction for multiple comparisons was performed than the prioritized gene set enrichment analysis using GENE2FUN from FUMA (Method 4F and Fig. 4). We identified that six gene set pathways had 18 gene set-PSC pairwise associations with 17 PSCs after Bonferroni correction for the number of gene sets and PSCs (N=16,768 and C from 32 to 1024, p-value threshold: 8.54 > -log10[p-value] > 7.03) (Fig. 3C, SI eFile 10). These gene sets imply critical biological and molecular pathways that might shape brain morphological changes and development. The reelin signaling pathway regulates neuronal migration, dendritic growth, branching, spine formation, synaptogenesis, and synaptic plasticity. [START_REF] Jossin | Reelin Functions, Mechanisms of Action and Signaling Pathways During Brain Development and Maturation[END_REF] The appendage morphogenesis and development pathways indicate how the anatomical structures of appendages are generated, organized, and progressed over time, often related to the cell adhesion pathway. These pathways elucidate how cells or tissues can be organized to create a complex structure like the human brain. [START_REF] Gilbert | Morphogenesis and Cell Adhesion[END_REF] In addition, the integral component of the cytoplasmic side of the endoplasmic reticulum membrane is thought to form a continuous network of tubules and cisternae extending throughout neuronal dendrites and axons. [START_REF] Wu | Contacts between the endoplasmic reticulum and other membranes in neurons[END_REF] The DSCAM (Down syndrome cell adhesion molecule) pathway likely functions as a cell surface receptor mediating axon pathfinding. Related proteins are involved in hemophilic intercellular interactions. [START_REF] Ly | DSCAM Is a Netrin Receptor that Collaborates with DCC in Mediating Turning Responses to Netrin-1[END_REF] Lastly, Nikolsky et al. [START_REF] Nikolsky | Genome-wide functional synergy between amplified and mutated genes in human breast cancer[END_REF] defined genes from the breast cancer 20Q11 amplicon pathway that were involved in the brain might indicate the brain metastasis of breast cancer, which is usually a late event with deleterious effects on the prognosis. [START_REF] Bos | Genes that mediate breast cancer metastasis to the brain[END_REF] In addition, previous findings [START_REF] Lanni | Cancer and Alzheimer's disease inverse relationship: an age-associated diverging derailment of shared pathways[END_REF][START_REF] Shafi | Inverse relationship between Alzheimer's disease and cancer, and other factors contributing to Alzheimer's disease: a systematic review[END_REF] revealed an inverse relationship between Alzheimer's disease and breast cancer, which might indicate a close genetic relationship between the disease and brain morphological changes mainly affecting the entorhinal cortex and hippocampus (PSC: C128_3 in Fig. 4).

Illustrations of genetic loci and pathways forming two patterns of structural covariance

To illustrate how underlying genetic underpinnings might form a specific PSC, we showcased two PSCs: C32_4 for the superior cerebellum and C128_3 for the hippocampus-entorhinal cortex. The two PSCs were highly heritable and polygenic in our GWAS using the entire UKBB data (Fig. 4, N=33,541). We used the FUMA 26 online platform to perform SNP2GENE for annotating the mapped genes and GENE2FUNC for prioritized gene set enrichment analyses (Method 4F). The superior cerebellum PSC was associated with genomic loci that can be mapped to 85 genes, which were enriched in many biological pathways, including psychiatric disorders, biological processes, molecular functions, and cellular components (e.g., apoptotic process, axon development, cellular morphogenesis, neurogenesis, and neuro differentiation).

For example, apoptosis -the regulated cell destruction -is a complicated process that is highly involved in the development and maturation of the human brain and neurodegenerative diseases. [START_REF] Yuan | Apoptosis in the nervous system[END_REF] Neurogenesis -new neuron formation -is crucial when an embryo develops and continues in specific brain regions throughout the lifespan. [START_REF] Steiner | A fresh look at adult neurogenesis[END_REF] All significant results of this prioritized gene set enrichment analysis are presented in SI eFile 11.

For the hippocampus-entorhinal cortex PSC, we mapped 45 genes enriched in gene sets defined from GWAS Catalog, including Alzheimer's disease and brain volume derived from hippocampal regions. The hippocampus and medial temporal lobe have been robust hallmarks of Alzheimer's disease. [START_REF] De Flores | Medial Temporal Lobe Networks in Alzheimer's Disease: Structural and Molecular Vulnerabilities[END_REF] In addition, these genes were enriched in the breast cancer 20Q11 amplicon pathway [START_REF] Nikolsky | Genome-wide functional synergy between amplified and mutated genes in human breast cancer[END_REF] and the pathway of metastatic breast cancer tumors [START_REF] Ginestier | Prognosis and gene expression profiling of 20q13-amplified breast cancers[END_REF] , which might indicate a specific distribution of brain metastases: the vulnerability of medial temporal lobe regions to breast cancer, [START_REF] Bos | Genes that mediate breast cancer metastasis to the brain[END_REF] or highlight an inverse association between Alzheimer's disease and breast cancer. [START_REF] Lanni | Cancer and Alzheimer's disease inverse relationship: an age-associated diverging derailment of shared pathways[END_REF] Lastly, the nuclear membrane encloses the cell's nucleus -the chromosomes reside inside -which is critical in cell formation activities related to gene expression and regulation. To further support the overlapping genetic underpinnings between this PSC and Alzheimer's disease, we calculated the genetic correlation (rg = -0.28; p-value=0.01) using GWAS summary statistics from the hippocampus-entorhinal cortex PSC (i.e., 33,541 people of European ancestry) and a previous independent study of Alzheimer's disease 31 (i.e., 63,926 people of European ancestry) using LDSC. [START_REF] Bulik-Sullivan | LD Score regression distinguishes confounding from polygenicity in genome-wide association studies[END_REF] All significant results of this prioritized gene set enrichment analysis are presented in SI eFile 12.

Multi-scale patterns of structural covariance derive disease-related imaging signatures

We used the multi-scale PSCs from a diverse population to derive imaging signatures that reflect brain development, aging, and the effects of several brain diseases. We investigate the added value of the multi-scale PSCs as building blocks of imaging signatures for several brain diseases and risk conditions using linear support vector machines (SVM) (Method 5). [START_REF] Davatzikos | Machine learning in neuroimaging: Progress and challenges[END_REF] The aim is to harness machine learning to drive a clinically interpretable metric for quantifying an individuallevel risk to each disease category. To this end, we define the signatures as SPARE-X (Spatial PAtterns for REcognition) indices, where X is the disease. For instance, SPARE-AD captures the degree of expression of an imaging signature of AD-related brain atrophy, which has been shown

to offer diagnostic and prognostic value in prior studies. [START_REF] Davatzikos | Longitudinal progression of Alzheimer's-like patterns of atrophy in normal older adults: the SPARE-AD index[END_REF] The most discriminative indices in our samples were SPARE-AD and SPARE-MCI (Fig. 5, SI eTable 4a and eFigure 5). C=1024 achieved the best performance for the single-scale analysis (e.g., AD vs. controls; balanced accuracy: 0.90±0.02; Cohen's d: 2.50). Multi-scale representations derived imaging signatures that showed the largest effect sizes to classify the patients from the controls (Fig. 5) (e.g., AD vs. controls; balanced accuracy: 0.92±0.02; Cohen's d: 2.61). PSCs obtained better classification performance than both AAL (e.g., AD vs. controls; balanced accuracy: 0.82±0.02; Cohen's d: 1.81) and voxel-wise regional volumetric maps (RAVENS) [START_REF] Davatzikos | Voxel-based morphometry using the RAVENS maps: methods and validation using simulated longitudinal atrophy[END_REF] (e.g., AD vs. controls; balanced accuracy: 0.85±0.02; Cohen's d: 2.04) (SI eTable 4a and eFigure 5). Our classification results were higher than previous baseline studies [START_REF] Wen | Convolutional neural networks for classification of Alzheimer's disease: Overview and reproducible evaluation[END_REF][START_REF] Samper-González | Reproducible evaluation of classification methods in Alzheimer's disease: Framework and application to MRI and PET data[END_REF] ,

which provided an open-source framework to objectively and reproducibly evaluate AD classification. Using the same cross-validation procedure and evaluation metric, they reported the highest balanced accuracy of 0.87±0.02 to classify AD from healthy controls. Notably, our experiments followed good practices, employed rigorous cross-validation procedures, and avoided critical methodological flaws, such as data leakage or double-dipping (refer to critical reviews on this topic elsewhere [START_REF] Wen | Convolutional neural networks for classification of Alzheimer's disease: Overview and reproducible evaluation[END_REF][START_REF] Kriegeskorte | Circular analysis in systems neuroscience: the dangers of double dipping[END_REF] ).

To test the robustness of these SPARE indices, we performed leave-one-site-out analyses for SPARE-AD using the combined 2003 PSCs from all scales (SI eTable 4b). Overall, holding the ADNI data out as independent test data resulted in a lower balanced accuracy (0.88±0.02) compared to the other cases for AIBL (0.95±0.02) and PENN data (0.95±0.02). The mean balanced accuracy (0.91±0.02) aligns with the nested cross-validated results using the full sample (Fig. 5).

BRIDGEPORT: bridging knowledge across patterns of structural covariance, genomics, and clinical phenotypes

We integrated our experimental results and the MuSIC atlas into the BRIDGEPORT online web portal. This online tool allows researchers to interactively browse the MuSIC atlas in 3D, query our experimental results via variants or PSCs, and download the GWAS summary statistics for further analyses. In addition, we allow users to search via conventional brain anatomical terms (e.g., the right thalamus proper) by automatically annotating traditional anatomic atlas ROIs, specifically from the MUSE atlas 39 (SI eTable 5), to MuSIC PSCs based on their degree of overlaps (SI eFigure 6). Open-source software dedicated to image processing, [START_REF] Doshi | MUSE: MUlti-atlas region Segmentation utilizing Ensembles of registration algorithms and parameters, and locally optimal atlas selection[END_REF] genetic quality check protocols, MuSIC generation with sopNMF, and machine learning 36 is also publicly 373 available (see Code Availability for details). 374

Discussion

The current study investigates patterns of structural covariance in the human brain at multiple scales from a large population of 50,699 people and, importantly, a very diverse cohort allowing us to capture patterns of structural covariance emanating from normal and abnormal brain development and aging, as well as from several brain diseases. Through extensive examination of the genetic architecture of these multi-scale PSCs, we confirmed genetic hits from previous T1-weighted MRI GWAS and, more importantly, identified 617 novel genomic loci and molecular and biological pathways that collectively influence brain morphological changes and development over the lifespan. Using a hypothesis-free, data-driven approach to first derive these PSCs using brain MRIs, we then uncovered their genetic underpinnings and further showed their potential as building blocks to predict various diseases. All experimental results and code are encapsulated and publicly available in BRIDGEPORT for dissemination: https://www.cbica.upenn.edu/bridgeport/, to enable various neuroscience studies to investigate these structural covariance patterns in diverse contexts. Together, the current study highlighted the adoption of machine learning methods in brain imaging genomics and deepened our understanding of the genetic architecture of the human brain.

Our findings reveal new insights into genetic underpinnings that influence structural covariance patterns in the human brain. Brain morphological development and changes are largely polygenic and heritable, and previous neuroimaging GWAS has not fully uncovered this genetic landscape. In contrast, genetic variants, as well as environmental, aging, and disease effects, exert pleiotropic effects in shaping morphological changes in different brain regions through specific biological pathways. The mechanisms underlying brain structural covariance are not yet fully understood. They may involve an interplay between common underlying genetic factors, shared susceptibility to aging, and various brain pathologies, which affect brain growth or degeneration in coordinated brain morphological changes. [START_REF] Alexander-Bloch | Imaging structural co-variance between human brain regions[END_REF] Our data-driven, multi-scale PSCs identify the hierarchical structure of the brain under the principle of structural covariance and are associated with genetic factors at different levels, including SNPs, genes, and gene set pathways.

These 617 novel genomic loci, as well as those previously identified, collectively shape brain morphological changes through many key biological and molecular pathways. These pathways are widely involved in reelin signaling, apoptotic processes, axonal development, cellular morphogenesis, neurogenesis, and neuro differentiation, [START_REF] Yuan | Apoptosis in the nervous system[END_REF][START_REF] Steiner | A fresh look at adult neurogenesis[END_REF] which may collectively influence the formation of structural covariance patterns in the brain. Strikingly, pathways involved in breast cancer shared overlapping genetic underpinnings evidenced in our MAGMA-based and prioritized (GENE2FUNC) gene set enrichment analyses (Fig. 3C and Fig. 4), which included specific pathways involved in breast cancer and metastatic breast cancer tumors. One previous study showed that common genes might mediate breast cancer metastasis to the brain, [START_REF] Bos | Genes that mediate breast cancer metastasis to the brain[END_REF] and a later study further corroborated that the metastatic spread of breast cancer to other organs (including the brain) accelerated during sleep in both mouse and human models. [START_REF] Diamantopoulou | The metastatic spread of breast cancer accelerates during sleep[END_REF] We further showcased that this brain metastasis of breast cancer might be associated with specific neuropathologic processes, which were captured by PSCs data driven by Alzheimer's diseaserelated neuropathology. For example, the hippocampus-entorhinal cortex PSC (C128_3, Fig. 4) connected the bilateral hippocampus and medial temporal lobe -the salient hallmark of Alzheimer's disease. Our gene set enrichment analysis results further support this claim: the genes were enriched in the gene sets of Alzheimer's disease and breast cancer (Fig. 4). Previous research [START_REF] Lanni | Cancer and Alzheimer's disease inverse relationship: an age-associated diverging derailment of shared pathways[END_REF][START_REF] Shafi | Inverse relationship between Alzheimer's disease and cancer, and other factors contributing to Alzheimer's disease: a systematic review[END_REF] also found an inverse association between Alzheimer's disease and breast cancer. In addition, PSCs from the cerebellum were the most genetically influenced brain regions, consistent with previous neuroimaging GWAS. [START_REF] Zhao | Genome-wide association analysis of 19,629 individuals identifies variants influencing regional brain volumes and refines their genetic co-architecture with cognitive and mental health traits[END_REF][START_REF] Elliott | Genome-wide association studies of brain imaging phenotypes in UK Biobank[END_REF] The cerebral cortex has been thought to largely contribute to the unique mental abilities of humans. However, the cerebellum may also be associated with a much more comprehensive range of complex cognitive functions and brain diseases than initially thought. [START_REF] Barton | Rapid Evolution of the Cerebellum in Humans and Other Great Apes[END_REF] Our results confirmed that many genetic substrates might support different molecular pathways, resulting in cerebellar functional organization, high-order functions, and dysfunctions in various brain disorders.

The current work demonstrates that appropriate machine learning analytics can be used to shed new light on brain imaging genetics. Previous neuroimaging GWAS leveraged multimodal imaging-derived phenotypes from conventional brain atlases 4,5 (e.g., the AAL atlas). In contrast, multi-scale PSCs are purely data-driven and likely to reflect the dynamics of underlying normal and pathological neurobiological processes giving rise to structural covariance. The diverse training sample from which the PSCs were derived, including healthy and diseased individuals of a wide age range, enriched the diversity of such neurobiological processes influencing the PSCs.

In addition, modeling structural covariance at multiple scales (i.e., multi-scale PSCs) indicated that disease effects could be robustly and complementarily identified across scales (Fig. 5), concordant with the paradigm of multi-scale brain modeling. [START_REF] Betzel | Multi-scale brain networks[END_REF] Imaging signatures of brain diseases, derived via supervised machine learning models, were consistently more distinctive when formed from multi-scale PSCs than single-scale PSCs. Multivariate learning techniques have gained significant prominence in neuroimaging and have recently attracted considerable attention in the domain of imaging genomics. These methods have proven valuable for analyzing complex and high-dimensional data, facilitating the exploration of relationships between imaging features and genetic factors. For instance, the MOSTest, a multivariate GWAS approach, preserves correlation structure among phenotypes via permutation on each SNP and derives a genotype vector for testing the association across all phenotypes [START_REF] Van Der Meer | Understanding the genetic determinants of the brain with MOSTest[END_REF] . A separate study by Soheili-Nezhad et al. demonstrated that genetic components obtained through PCA or ICA applied to neuroimaging GWAS summary statistics exhibited greater reproducibility than raw univariate GWAS effect sizes [START_REF] Soheili-Nezhad | Reproducibility of Principal and Independent Genomic Components of Brain Structure and Function[END_REF] . A recent study utilized a CNN-based autoencoder to discover new phenotypes and identify numerous novel genetic signals [START_REF] Patel | New phenotype discovery method by unsupervised deep representation learning empowers genetic association studies of brain imaging[END_REF] . Despite the effectiveness of these multivariate approaches in GWAS, they typically conduct phenotype engineering before performing GWAS without explicitly incorporating imaging genetic associations during the modeling process. Yang et al. recently conducted a study that employed generative adversarial networks (termed GeneSGAN 45 ) to integrate imaging and genetic variations within the modeling framework to address this limitation. By incorporating both modalities, their approach aimed to capture the complexity and heterogeneity of disease manifestations.

MuSIC -with the strengths of being data-driven, multi-scale, and disease-effect informative -contributes to the century-old quest for a "universal" atlas in brain cartography [START_REF] Eickhoff | Imaging-based parcellations of the human brain[END_REF] and is highly complementary to previously proposed brain atlases. For instance, Chen and colleagues 47 used a semi-automated fuzzy clustering technique with MRI data from 406 twins and parcellated the cortical surface area into a genetic covariance-informative brain atlas; MuSIC was data-driven by structural covariance. Glasser and colleagues 48 adopted a semi-automated parcellation procedure to create a multimodal cortex atlas from 210 healthy individuals.

Although this method successfully integrates multimodal information from cortical folding, myelination, and functional connectivity, this semi-automatic approach requires significant resources, some with limited resolution. MuSIC allows flexible, multiple scales for delineating macroscopic brain topology; including patient samples exposes the model to sources of variability that may not be visible in healthy controls. Another pioneering endeavor is the Allen Brain Atlas project, [START_REF] Sunkin | Brain Atlas: an integrated spatio-temporal portal for exploring the central nervous system[END_REF] whose overarching goals of mapping the human brain to gene expression data via existing conventional atlases, identifying local gene expression patterns across the brain in a few individuals, and deepening our understanding of the human brain's differential genetic architecture, are complementary to ours -characterizing the global genetic architecture of the human brain, emphasizing pathogenic variability and morphological heterogeneity.

Bridging knowledge across the brain imaging, genomics, and machine learning communities is another pivotal contribution of this work. BRIDGEPORT provides a platform to lower the entry barrier for whole-brain genetic-structural analyses, foster interdisciplinary communication, and advocate for research reproducibility. [START_REF] Wen | Convolutional neural networks for classification of Alzheimer's disease: Overview and reproducible evaluation[END_REF][START_REF] Munafò | A manifesto for reproducible science[END_REF][START_REF] Poldrack | Scanning the horizon: towards transparent and reproducible neuroimaging research[END_REF][START_REF] Routier | Clinica: An Open-Source Software Platform for Reproducible Clinical Neuroscience Studies[END_REF][START_REF] Wen | Reproducible Evaluation of Diffusion MRI Features for Automatic Classification of Patients with Alzheimer's Disease[END_REF] The current study demonstrates the broad applicability of this large-scale, multi-omics platform across a spectrum of neurodegenerative and neuropsychiatric diseases.

The present study has certain limitations. Firstly, the sopNMF method utilized in brain parcellation considers only imaging structural covariance and overlooks the genetic determinants contributing to forming these structural networks, as indicated by our GWAS findings.

Consequently, further investigations are needed to integrate imaging and genetics into brain parcellation. Additionally, it is important to note that our GWAS analyses primarily involved participants of European ancestry. To enhance genetic findings for underrepresented ethnic groups, future studies should prioritize the inclusion of diverse ancestral backgrounds, thereby promoting a more comprehensive understanding of the genetic underpinnings across different populations.

Methods

Method 1: Structural covariance patterns via stochastic orthogonally projective nonnegative matrix factorization

The sopNMF algorithm is a stochastic approximation built and extended based on opNMF [START_REF] Sotiras | Finding imaging patterns of structural covariance via Non-Negative Matrix Factorization[END_REF][START_REF] Yang & Oja | Linear and Nonlinear Projective Nonnegative Matrix Factorization[END_REF] .

We consider a dataset of 𝑛 MR images and 𝑑 voxels per image. We represent the data as a matrix X where each column corresponds to a flattened image:

𝑿 = [𝑥 ! , 𝑥 " , … , 𝑥 # ], 𝑿 ∈ ℝ $% &×# .
The

𝑿 𝒃 ∈ ℝ $% &×-of 𝑏 ≤
𝑛 images (b is the batch size; b=32 was used in the current analyses); this is done without replacement so that all data goes through the model once (⌈𝑛/𝑏⌉). In this case, the updating rule can be rewritten as

𝑾 ./! = 𝑾 . O𝑿 𝒃 𝑿 𝒃 𝑻 𝑾P . (𝑾𝑾 0 𝑿 𝒃 𝑿 𝒃 𝑻 𝑾) . ( 2
)
We calculate the loss on the entire dataset at the end of each epoch (i.e., the loss is incremental across all batches) with the following expression:

R S𝑿 𝒃_𝒊 -𝑾𝑾 𝑻 𝑿 𝒃𝒊 S ) " ⌈#/-⌉ +6! (3)
We evaluated the training loss and the sparsity of W at the end of each iteration. Moreover, early stopping was implemented to improve training efficiency and alleviate overfitting. We summarize the sopNMF algorithm in SI Algorithm 1. An empirical comparison between sopNMF and opNMF is detailed in SI eMethod 1.

We applied sopNMF to the training population (N=4000). The component matrix W was sparse after the algorithm converged with a pre-defined maximum number of epochs (100 by default) with an early stopping criterion. To build the MuSIC atlas, we clustered each voxel (row-wise) into one of the 𝑟 features/PSCs as follows:

𝑴 7 = argmax 8 (𝑾 7,8 ) (4)
where M is a d-dimensional vector and 𝑗 ∈ {1. . 𝑑}. The j-th element of M equals k if 𝑾 7,8 is the maximum value of the j-th row. Intuitively, M indicates which of the r PSCs each voxel belongs to. We finally projected the vector 𝑴 ∈ ℝ $% & into the original image space to visualize each PSC of the MuSIC atlas (Fig. 1). Of note, 13 PSCs have vanished in this process for C=1024: all 0 for these 13 vectors.

Method 2: Study population

We consolidated a large-scale multimodal consortium (N=50,699) consisting of imaging, cognition, and genetic data from 12 studies, 130 sites, and 12 countries. We present the detailed demographic information of the population under study in SI eTable 1. All individual studies were approved by their local corresponding Institutional Review Boards (IRB) (SI eText 2). This large-scale consortium reflects the diversity of MRI scans over different races, disease conditions, and ages over the lifespan. To be concise, we defined four populations or data sets per analysis across the paper: i) discovery set, ii) replication set, iii) training population, and iv) comparison population (refer to SI eText 3 for details).

Method 3: Image processing and statistical harmonization

(A): Image processing. Images that passed the quality check (SI eMethod 4) were first corrected for magnetic field intensity inhomogeneity. [START_REF] Tustison | N4ITK: improved N3 bias correction[END_REF] Voxel-wise regional volumetric maps (RAVENS) [START_REF] Davatzikos | Voxel-based morphometry using the RAVENS maps: methods and validation using simulated longitudinal atrophy[END_REF] for each tissue volume were then generated by using a registration method to spatially align the skull-stripped images to a template in MNI-space. [START_REF] Ou | Deformable Registration via Attribute Matching and Mutual-Saliency Weighting[END_REF] We applied sopNMF to the RAVENS maps to derive MuSIC.

(B): Statistical harmonization of MuSIC PSCs:

We applied MuSIC to the entire population (N=50,699) to extract the multi-scale PSCs. Specifically, MuSIC was applied to each individual's RAVENS gray matter map to extract the sum of brain volume in each PSC. Subsequently, the PSCs were statistically harmonized by an extensively validated approach, i.e., ComBat-GAM [START_REF] Pomponio | Harmonization of large MRI datasets for the analysis of brain imaging patterns throughout the lifespan[END_REF] (SI eMethod 3) to account for site-related differences in the imaging data. After harmonization, the PSCs were normally distributed (skewness = 0.11±0.17, and kurtosis = 0.67±0.68) (SI eFigure 7A and B). To alleviate the potential violation of normal distribution in downstream statistical learning, we quantile-transformed all PSCs. In agreement with the literature, 57,58 males were found to have larger brain volumes than females on average (SI eFigure 7C). Overall, the (B): Heritability estimates and genome-wide association analysis. We estimated the SNPbased heritability explained by all autosomal genetic variants using GCTA-GREML. [START_REF] Yang | GCTA-GREML accounts for linkage disequilibrium when estimating genetic variance from genome-wide SNPs[END_REF] We adjusted for confounders of age (at imaging), age-squared, sex, age-sex interaction, age-squaredsex interaction, ICV, and the first 40 genetic principal components (PC), guided by a previous neuroimaging GWAS [START_REF] Zhao | Genome-wide association analysis of 19,629 individuals identifies variants influencing regional brain volumes and refines their genetic co-architecture with cognitive and mental health traits[END_REF] . In addition, Elliot et al. [START_REF] Elliott | Genome-wide association studies of brain imaging phenotypes in UK Biobank[END_REF] investigated more than 200 confounders in another study. Therefore, our sensitivity analyses included four additional imaging-related covariates (i.e., brain positions and head motion). One-side likelihood ratio tests were performed to derive the heritability estimates. In GWAS, we performed a linear regression for each PSC and included the same covariates as in the heritability estimates using PLINK. for the definition of the index, candidate, independent significant, lead SNP, and genomic locus).

In particular, the threshold for significance was set to 5×10 -8 (clump-p1) for the index SNPs and 0.05 (clump-p2) for the candidate SNPs. The threshold for linkage disequilibrium-based clumping was set to 0.60 (clump-r2) for independent significant SNPs and 0.10 for lead SNPs.

The linkage disequilibrium physical-distance threshold was 250 kilobases (clump-kb). Genomic loci consider linkage disequilibrium (within 250 kilobases) when interpreting the association results. The GWASRAPIDD 65 package (version: 0.99.14) was then used to query the genomic loci for any previously-reported associations with clinical phenotypes documented in the NHGRI-EBI GWAS Catalog 15 (p-value < 1.0×10 -5 , default inclusion value of GWAS Catalog).

significant chromatin interaction between the disease-associated regions and nearby or distant genes. [START_REF] Watanabe | Functional mapping and annotation of genetic associations with FUMA[END_REF] In addition, GENE2FUNC studies the expression of prioritized genes and tests for the enrichment of the set of genes in pre-defined pathways. We used the mapped genes as prioritized genes. The background genes were specified as all genes in FUMA, and all other parameters were set by default. We only reported gene sets with adjusted p-value < 0.05.

(G): Correction for multiple comparisons. We practiced a conservative procedure to control for the multiple comparisons. In the case of GWAS, we chose the default genome-wide significant threshold (5.0x10 -8 , and 0.05 for all other analyses) and independently adjusted for multiple comparisons (Bonferroni methods) at each scale by the number of PSCs. We corrected the p-values for the number of phenotypes (N=6) for genetic correlation analyses. We adjusted the p-values for the number of PSCs at each scale for heritability estimates. For gene analyses, we controlled for both the number of PSCs at each scale and the number of genes. We adopted these strategies per analysis to correct the multiple comparisons because PSCs of different scales are likely hierarchical and correlated -avoiding the potential of "overcorrection".

(H): Replication analysis for genome-wide association studies. We performed GWAS by fitting the same linear regressing models as the discovery set. Also, following the same procedure for consistency, we corrected the multiple comparisons using the Bonferroni method.

We corrected it for the number of genomic loci (N=915) found in the discovery set with a nominal p-value of 0.05, which thereby resulted in a stringent test with an equivalent p-value threshold of 3.1x10 -5 (i.e., (-log10[p-value] = 4.27). We performed a replication for the 915 genomic loci, but, in reality, SNPs in linkage disequilibrium with the genomic loci are likely highly significant. We demonstrate how underlying genomic loci and biological pathways might influence the formation, development, and changes of two specific PSCs: the 4 th PSC of the C32 PSCs (C32_4) that resides in the superior part of the cerebellum and the 3 rd PSC of the C128 PSCs (C128_3) that includes the bilateral hippocampus and entorhinal cortex. We first performed SNP2GENE to annotate the mapped genes in the Manhattan plots and then ran GENE2FUNC for the prioritized gene set enrichment analysis (Method 4F). The mapped genes are input genes for prioritized gene set enrichment analyses. The heat map shows the significant gene sets from the GWAS Catalog, curated genes, and gene ontology (GO) that survived the correction for multiple comparisons. We selectively present the schematics for three pathways: apoptosis, neurogenesis, and nuclear membrane function. Several other key pathways are highlighted in bold, and the 3D maps of the two PSCs are presented. Imaging signatures (SPARE indices) of brain diseases, derived via supervised machine learning models, are more distinctive when formed from multi-scale PSCs than single-scale PSCs. The kernel density estimate plot depicts the distribution of the patient group (blue) in comparison to the healthy control group (red), reflecting the discriminative power of the diagnosis-specific SPARE (imaging signature) indices. We computed Cohen's d for each SPARE index between groups to present the effect size of its discrimination power. * represents the model with the largest Cohen's d for each SPARE index to separate the control vs. patient groups; # represents the model with the best performance with single-scale PSCs. Our results demonstrate that the multi-scale PSCs generally achieve the largest discriminative effect sizes (ES) (SI eTable 4a). As a reference, Cohen's d of ≥ 0.2, ≥ 0.5, and ≥ 0.8, respectively, refer to small, moderate, and large effect sizes.

  sopNMF algorithm factorizes X into two low-rank (𝑟) matrices 𝑾 ∈ ℝ $% &×( and 𝑯 ∈ ℝ $% (×# under the constraints of non-negativity and column-orthonormality. Using the Frobenius norm, the loss of this factorization problem can be formulated as ‖𝑿 -𝑾𝑯‖ ) " subject to 𝑯 = 𝑾 𝑻 𝑿 , 𝑾 ≥ 0 and 𝑾 𝑻 𝑾 = 𝑰 (1)where I stands for the identity matrix. The columns 𝑤 + ∈ ℝ & , ‖𝑤 + ‖ " = 1, ∀ 𝑖 ∈ {1. . 𝑟} of the socalled component matrix 𝑾 = [𝑤 ! , 𝑤 " , … , 𝑤 ( ] are part-based representations promoting sparsity in data in this lower-dimensional subspace. From this perspective, the loading coefficient matrix 𝑯 represents the importance (weights) of each feature above for a given image. Instead of optimizing the non-convex problem in a batch learning paradigm (i.e., reading all images into memory) as opNMF, 9 sopNMF subsamples the number of images at each iteration, thereby significantly reducing its memory demand by randomly drawing data batches

64

 64

(

  C): Identification of novel genomic loci. Using PLINK, we clumped the GWAS summary statistics based on their linkage disequilibrium to identify the genomic loci (see SI eMethod 5

Method 5 : 2 .

 52 Pattern analysis via machine learning for individualized imaging signatures SPARE-AD captures the degree of expression of an imaging signature of AD, and prior studies have shown its diagnostic and prognostic values.[START_REF] Davatzikos | Longitudinal progression of Alzheimer's-like patterns of atrophy in normal older adults: the SPARE-AD index[END_REF] Here, we extended the concept of the SPARE imaging signature to multiple diseases (SPARE-X, X represents disease diagnoses). Following our reproducible open-source framework[START_REF] Samper-González | Reproducible evaluation of classification methods in Alzheimer's disease: Framework and application to MRI and PET data[END_REF] , we performed nested cross-validation (SI eMethod 6) for the machine learning models and derived imaging signatures to quantify individualized disease vulnerability.SPARE indices.MuSIC PSCs were fit into a linear support vector machine (SVM) to derive SPARE-AD, MCI, SCZ, DM, HTN, MDD, and ASD. Specifically, the SVM aims to classify the patient group (e.g., AD) from the control group and outputs a continuous variable (i.e., the SPARE indices), which indicates the proximity of each participant to the hyperplane in either the patient or control space. We compared the classification performance using different sets of features: i) the single-scale PSC from 32 to 1024, ii) the multi-scale PSCs by combining all features (with and without feature selections embedded in the CV); iii) the ROIs from the AAL atlas; and iv) voxel-wise RAVENS maps. The samples selected for each task are presented in SI eTable No statistical methods were used to predetermine the sample size. The experiments were not randomized, and the investigators were not blinded to allocation during experiments and outcome assessment.

Figures Figure 1 :AC

 1 Figures

Figure 3 :

 3 Figure 3: Patterns of structural covariance highlight novel genomic loci and pathways that shape the human brain. A) Patterns of structural covariance (PSC) in the human brain are polygenic: the number of genomic loci of each PSC is projected onto the image space to show a statistical brain map characterized by the number (C) of PSCs. In addition, common genetic variants exert pleiotropic effects on the PSCs: circular plots showed the number of associated PSCs (histograms in blue color) of each genomic loci over the entire autosomal chromosome (1-22). The histogram was plotted for the number of PSCs for each genomic locus in the circular plots. B) Novel genomic loci revealed by the multi-scale PSCs compared to previous findings from the GWAS Catalog, 15 T1-weighted MRI GWAS 4,5 , and the AAL atlas regions of interest. The green bar indicates the 617 novel genomic loci not previously associated with any clinical traits in GWAS Catalog; the black bar presents the loci identified in other studies that overlap (grey bar for loci in linkage disequilibrium) with the loci from our results; the yellow bar indicates the unique loci in other studies. C) Pathway enrichment analysis highlights six unique biological pathways and functional categories (after Bonferroni correction for 16,768 gene sets and the number of PSCs) that might influence the changes of PSCs. DSCAM: Down syndrome cell adhesion molecule.

2 Figure 4 :

 24 Figure 4: Illustrations of multiple genetic loci and pathways shaping specific patterns of structural covarianceWe demonstrate how underlying genomic loci and biological pathways might influence the formation, development, and changes of two specific PSCs: the 4 th PSC of the C32 PSCs (C32_4) that resides in the superior part of the cerebellum and the 3 rd PSC of the C128 PSCs (C128_3) that includes the bilateral hippocampus and entorhinal cortex. We first performed SNP2GENE to annotate the mapped genes in the Manhattan plots and then ran GENE2FUNC for the prioritized gene set enrichment analysis (Method 4F). The mapped genes are input genes for prioritized gene set enrichment analyses. The heat map shows the significant gene sets from the GWAS Catalog, curated genes, and gene ontology (GO) that survived the correction for multiple comparisons. We selectively present the schematics for three pathways: apoptosis, neurogenesis, and nuclear membrane function. Several other key pathways are highlighted in bold, and the 3D maps of the two PSCs are presented.

47 Figure 5 :

 475 Figure5: Individualized imaging signatures based on pattern analysis via machine learning. Imaging signatures (SPARE indices) of brain diseases, derived via supervised machine learning models, are more distinctive when formed from multi-scale PSCs than single-scale PSCs. The kernel density estimate plot depicts the distribution of the patient group (blue) in comparison to the healthy control group (red), reflecting the discriminative power of the diagnosis-specific SPARE (imaging signature) indices. We computed Cohen's d for each SPARE index between groups to present the effect size of its discrimination power. * represents the model with the largest Cohen's d for each SPARE index to separate the control vs. patient groups; # represents the model with the best performance with single-scale PSCs. Our results demonstrate that the multi-scale PSCs generally achieve the largest discriminative effect sizes (ES) (SI eTable 4a). As a reference, Cohen's d of ≥ 0.2, ≥ 0.5, and ≥ 0.8, respectively, refer to small, moderate, and large effect sizes.

Table 1

 1 provides an overview of the abbreviations used in the present study.

Table 1 . Abbreviations used in the present study 190

 1

	Item	Abbreviation	Item	Abbreviation
	Pattern of structural covariation	PSC	Independent component analysis	ICA
	Genome-wide association study	GWAS	BRaIn knowleDGE PORTal	BRIDGEPORT
	Orthogonal projective non-negative matrix factorization	opNMF	Multi-scale Structural Imaging Covariance	MuSIC
	Stochastic orthogonal
	projective non-negative	sopNMF	Machine learning	ML
	matrix factorization
	Principal component analysis	PCA	UK Biobank	UKBB
	Imaging-based		Psychosis
	coordinate SysTem for		Heterogeneity
	AGing and	iSTAGING	Evaluated via	PHENOM
	NeurodeGenerative		Dimensional
	diseases		Neuroimaging
	Single nucleotide polymorphism	SNP	Region of interest	ROI
	Magnetic resonance imaging	MRI	Automated anatomical labeling	AAL
	MUlti-atlas region
	Segmentation utilizing	MUSE	Alzheimer's disease	AD
	Ensembles
	Spatial PAtterns for REcognition	SPARE	Support vector machine	SVM

Method 3B). Unit II (Fig. 1B) investigates the harmonized data for 2003 PSCs (13 PSCs have vanished in this process for C=1024; see Method 1) in two brain structural covariance analyses. Specifically, we performed i) GWAS (Method 4)

  [START_REF] Seeley | Neurodegenerative diseases target large-scale human brain networks[END_REF] Multi-scale PSCs were extracted across the entire population and statistically harmonized[START_REF] Pomponio | Harmonization of large MRI datasets for the analysis of brain imaging patterns throughout the lifespan[END_REF] (

617 novel genomic loci of patterns of structural covariance

	sets) in the UK Biobank (UKBB) data. The scatter plot of the two sets of h 2 estimates is shown in
	SI eFigure 3. The h 2 estimates and p-values for all PSCs are detailed in SI eFile 1 (discovery set)
	and eFile 2 (replication set). Our results confirm that brain structure is heritable to a large extent
	and identify the spatial distribution of the most highly heritable regions of the brain (e.g.,
	subcortical gray matter structures and cerebellum regions). 14
	We discovered genomic locus-PSC pairwise associations (Method 4C, SI eMethod 5) within
	the discovery set and then independently replicated these associations on the replication set. We
	found that 915 genomic loci had 3791 loci-PSC pairwise significant associations with 924 PSCs
	after Bonferroni correction (Method 4G) for the number of PSCs (p-value threshold per scale:
	10.3 > -log10[p-value] > 8.8) (SI eFile 3, and Fig. 3A). Our results showed that the formation of
	these PSCs is largely polygenic; the associated SNPs might play a pleiotropic role in shaping
	these networks.
	Compared to previous literature, out of the 915 genomic loci, the multi-scale PSCs
	identified 617 novel
	14, 0.29±0.15, 0.25±0.15, 0.27±0.15, 0.31±0.15 for scales C=32, 64, 128,
	256, 512 and 1024 of the PSCs, respectively. The Pearson correlation coefficient between the two
	independent estimates of h 2 was r = 0.94 (p-value < 10 -6 , between the discovery and replication

genomic loci not previously associated with any traits or phenotypes in the GWAS Catalog 15 (SI eFile 4, Fig. 3B, query date: April 5 th

	, 2023). These novel associations
	might indicate subtle neurobiological processes that are captured thanks to the biologically
	relevant structural covariance expressed by sopNMF. The multi-scale PSCs identified many
	novel associations by constraining this comparison to previous neuroimaging GWAS 12,13 using
	T1w MRI-derived phenotypes (e.g., regions of interest from conventional brain atlases) (Fig 3B,

SI eTable 3, eFile 5, 6, and 7). Our

UKBB replication set analysis (Method 4H) demonstrated that 3638 (96%) exact genomic locus-PSC associations were replicated at nominal significance (-log10[p-value] > 1.31), 2705 (72%) of which were significant after correction for multiple comparisons (Method 4G, -log10[p-value] > 4.27). We present this validation in SI eFile 8 from the replication set. The summary statistics, Manhattan, and QQ plots derived from the combined population (N=33,541)

Acknowledgments

We acknowledge the contribution from the iSTAGING, the BLSA, the BIOCARD, the PHENOM, the ADNI studies, and the AI4AD consortium. The initial funding package for WJ as an Assistant Professor of Neurology, provided by Stevens Neuroimaging and Informatics Institute, Keck School of Medicine of USC, University of Southern California, supports the

Data Availability

The GWAS summary statistics corresponding to this study are publicly available on the BRIDGEPORT web portal (https://www.cbica.upenn.edu/bridgeport/) and the MEDICINE web portal (http://labs.loni.usc.edu/medicine/).

Code Availability

The software and resources used in this study are all publicly available:

• sopNMF: https://pypi.org/project/sopnmf/, MuSIC, and sopNMF (developed for this study)

• BRIDGEPORT: https://www.cbica.upenn.edu/bridgeport/, (developed for this study)

• MLNI: https://pypi.org/project/mlni/, machine learning (developed for this study) • MUSE: https://www.med.upenn.edu/sbia/muse.html, image preprocessing • PLINK: https://www.cog-genomics.org/plink/, GWAS • GCTA: https://yanglab.westlake.edu.cn/software/gcta/#Overview, heritability estimates • LDSC: https://github.com/bulik/ldsc, genetic correlation estimates • MAGMA: https://ctg.cncr.nl/software/magma, gene analysis • GWASRAPIDD: https://rmagno.eu/gwasrapidd/articles/gwasrapidd.html, GWAS Catalog query • MsigDB: https://www.gsea-msigdb.org/gsea/msigdb/, gene sets database

Combat-GAM model slightly improved data normality across sites (SI eFigure 7E-H). The AAL ROIs underwent the same statistical harmonization procedure.

Method 4: Genetic analyses

Genetic analyses were restricted to the discovery and replication set from UKBB (Method 2).

We processed the array genotyping and imputed genetic data (SNPs). The two data sets went through a "best-practice" imaging-genetics quality check (QC) protocol (Method 4A) and were restricted to participants of European ancestry. This resulted in 18,052 participants and 8,430,655 SNPs for the discovery set and 15,243 participants and 8,470,709 SNPs for the replication set.

We reperformed the genetic QC and genetic analyses for the combined populations for BRIDGEPORT, resulting in 33,541 participants and 8,469,833 SNPs. Method 4G details the correction for multiple comparisons throughout our analyses.

(A): Genetic data quality check protocol. First, we excluded related individuals (up to 2 nddegree) from the complete UKBB sample (N=488,377) using the KING software for family relationship inference. [START_REF] Manichaikul | Robust relationship inference in genome-wide association studies[END_REF] We then removed duplicated variants from all 22 autosomal chromosomes. We also excluded individuals for whom either imaging or genetic data were not available. Individuals whose genetically identified sex did not match their self-acknowledged sex were removed. Other excluding criteria were: i) individuals with more than 3% of missing genotypes; ii) variants with minor allele frequency (MAF) of less than 1%; iii) variants with larger than 3% missing genotyping rate; iv) variants that failed the Hardy-Weinberg test at 1x10 -10 . To adjust for population stratification, [START_REF] Price | New approaches to population stratification in genome-wide association studies[END_REF] we derived the first 40 genetic principle components (PC) using the FlashPCA software [START_REF] Abraham | FlashPCA2: principal component analysis of Biobankscale genotype datasets[END_REF] . The genetic pipeline was also described elsewhere [START_REF] Wen | Characterizing Heterogeneity in Neuroimaging, Cognition, Clinical Symptoms, and Genetics Among Patients With Late-Life Depression[END_REF] .

We defined a genomic locus as novel when it was not present in GWAS Catalog (query date: April 5 th , 2023).

(D): Gene-level associations with MAGMA. We performed gene-level association analysis using MAGMA. [START_REF] Leeuw | MAGMA: Generalized Gene-Set Analysis of GWAS Data[END_REF] (E): Hypothesis-free gene set enrichment analysis with MAGMA. Using the gene-level association p-values, we performed gene set enrichment analysis using MAGMA. Gene sets were obtained from Molecular Signatures Database (MsigDB, v7.5.1), [START_REF] Subramanian | Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles[END_REF] including 6366 curated gene sets and 10,402 Gene Ontology (GO) terms. All other parameters were set by default for MAGMA. This hypothesis-free analysis resulted in a more stringent correction for multiple comparisons (i.e., by the total number of tested genes and PSCs) than the FUMA-prioritized gene set enrichment analysis (see below F).

(F): FUMA analyses for the illustrations of specific PSCs. In SNP2GENE, three different methods were used to map the SNPs to genes. First, positional mapping maps SNPs to genes if the SNPs are physically located inside a gene (a 10 kb window by default). Second, expression quantitative trait loci (eQTL) mapping maps SNPs to genes showing a significant eQTL association. Lastly, chromatin interaction mapping maps SNPs to genes when there is a MZ serves as a consultant and/or speaker for the following pharmaceutical companies: Eurofarma, Lundbeck, Abbott, Greencare, Myralis, and Elleven Healthcare.

Competing Interests

Authors' contributions

Dr. Wen takes full responsibility for the integrity of the data and the accuracy of the data analysis.