A Network Approach to Mapping Mouse Brain-wide Mitochondrial Respiratory Chain Capacity in Relation to Behavior

1 Division of Behavioral Medicine, Department of Psychiatry, Columbia University Irving Medical Center, New York NY 2 Department of Psychiatry and Behavioral Sciences, Stanford University 3 The Sackler Institute, Department of Psychiatry, Columbia University Irving Medical Center, New York NY 4 Department of Biological Sciences, Columbia University 5 Division of Systems Neuroscience, Department of Psychiatry, Columbia University Irving Medical Center, New York NY 6 Brain Mind Institute, Ecole Polytechnique Federal de Lausanne (EPFL), Switzerland Division of Molecular Therapeutics, Department of Psychiatry, Columbia University Irving Medical Center, New York NY 8 Department of Pediatrics, Columbia University Irving Medical Center, New York NY 9 Division of Developmental Neuroscience, Department of Psychiatry, Columbia University Irving Medical Center, New York NY 10 New York State Psychiatric Institute, New York NY 11 Department of Neurology, Columbia University Irving Medical Center, New York NY


Introduction
The shaping of behaviors by life experiences is driven by energetically demanding circuitry across the brain 1 . The brain's enormous energetic demand is mainly met by ATP produced through oxidative phosphorylation (OxPhos), subserved by the combined activities of respiratory chain (RC) enzymes within mitochondria 2 . As a result, mitochondria influence multiple aspects of brain development and function ranging from dendritic and axonal branching 3,4 , remodeling of gene expression 5 , regulation of neurotransmitter release and excitability in mature synapses [6][7][8] , neurogenesis 9 , and inflammation 10 . Thus, the mounting molecular and functional evidence that the brain is under energetic constraints 11,12 suggests that if we want to understand the basis of brain function and behavior, we must understand key aspects of brain mitochondrial biology.
Mitochondria are small, dynamic, multifunctional organelles with their own genome 13 , but not all mitochondria are created equal. Mitochondria serving different cellular demands (i.e., in different cell types) have different relative molecular compositions, morphologies, and functional phenotypes [14][15][16][17][18][19] . Therefore, developing a comprehensive understanding of the association between mitochondrial biology and animal behavior calls for assessments of multiple functional and molecular mitochondrial features, across multiple brain regions. By analyzing mitochondrial features across multiple brain regions simultaneously, we can also potentially uncover unknown brain "mitochondrial circuits". This idea aligns with the evolving understanding of large-scale brain circuitry and metabolism 11,20 , and of the network distribution of neural activity achieved using brain-wide, high spatial resolution methods (e.g., MRI-based functional and structural connectivity maps) 21,22 . Although a similar degree of resolution for mitochondrial phenotyping is not feasible with current technologies, miniaturization of biochemical and molecular mitochondrial assays open new possibilities to systematically map mitochondria-to-behavior associations across multiple cortical and sub-cortical brain regions in the same animal.
Both mitochondrial biology and animal behaviors measured on standardized tests exhibit naturally occurring and acquired (i.e., experience-dependent) variation, providing an opportunity to map their associations. For example, behaviorally, exposure to social stress such as chronic social defeat predictably increases potential threat vigilance and social avoidance 23 , but animals differ in their vulnerability/resilience to these stress-induced behavior changes 24,25 . Experimental challenges aimed at modeling neuroendocrine disturbances resulting from chronic stress, such as chronic exposure to corticosterone, also induce avoidance behaviors associated with recalibrations in specific brain circuitry 26 , gene expression 27,28 , and anatomical plasticity (i.e., atrophy) in stress-sensitive brain regions like the hippocampus 29 . In relation to mitochondria, a separate body of literature similarly documents naturally-occurring mitochondrial variation [30][31][32] , as well as stress-induced functional mitochondrial recalibrations that occur within days to weeks of stress exposure [33][34][35] (meta-analysis in 36 ), potentially linking mitochondrial biology to behavior.
More direct causal experiments show that mitochondrial RC enzyme activities directly influence the brain and specific behavioral domains including working memory 37,38 , social dominance 31,39,40 and anxiety-related behavior 4 . Targeted mitochondrial defects even cause mood disorder-like phenotypes in animals 41 , positioning mitochondria as upstream modulators of brain function and behavior. Moreover, mitochondrial RC defects are likely implicated in the etiology of psychiatric, neurological and degenerative disorders in humans 42 , and in vivo brain metabolic imaging studies show that energy metabolism in specific brain areas (e.g., nucleus accumbens, NAc) predict cognitive performance and anxiety [43][44][45] , making these biological questions also potentially relevant to human mental health.
Although the importance of mitochondria for brain structure and function is unequivocal, we lack an understanding of potential differences in mitochondria across different brain regions.
This calls for more systematic mapping of region-specific brain mitochondrial biology in relation to behavior. To address the hypothesis that mitochondrial phenotypes in specific brain regions are associated with behaviors, we have further miniaturized existing biochemical and molecular assays of mitochondrial OxPhos enzyme activities for sub-milligram tissue samples and deployed them across 17 cortical and sub-cortical brain regions in mice with a wide range of behavioral phenotypes, quantified through four behavioral tests. Using network-based connectivity analysis, we also sought to explore the distribution of mitochondrial phenotypes across brain regions, finding evidence that mouse brain mitochondria specialize as distinct functional networks linked to behavior. Further, we use gene co-expression and structural/anatomical connectivity data, which anchor the newly observed mitochondrial circuits into other modalities, and provide a foundation for future mechanistic studies to elucidate the specific energetic contributions to behavioral variation.

Miniaturization of mitochondrial assays for sub-milligram resolution
Although enzymatic activity assays directly reflect the capacity of RC complexes and therefore energy production capacity, previously available assays suffer from low throughput that preclude a multi-region, brain-wide analysis in dozens of animals. Building from our efforts to miniaturize and scale the throughput of mitochondrial RC activity assays in immune cells 30,32 , here we miniaturized and optimized spectrophotometric assays in 96-well plate format for brain tissue, validated against the standard cuvette-based reaction (Extended Data Fig. 1, and Methods). Using this optimized platform, we can quantify the enzymatic activities of RC complex I (CI, NADH-ubiquinone oxidoreductase), complex II (CII, succinate-ubiquinone oxidoreductase), complex IV (CIV, cytochrome c oxidase), and citrate synthase (CS, a Krebs cycle enzyme and marker of mitochondrial content) in two 1mm-diameter, 200μm deep tissue punches, establishing in our hands the lowest detection limit for mouse brain tissue. This represents <1mg of tissue (estimated 0.33mg), over an order of magnitude more sensitive than currently available methods.
In this pipeline, the same biological sample used for enzymatic activities is also used to quantify mtDNA abundance, both by i) the classical metric relative to the nuclear genome -mtDNA/nDNA ratio, termed mtDNA copy number (mtDNAcn) 46 , and ii) per tissue volume (mtDNA copies per μ m 3 or per mg), termed mtDNA density. qPCR data of nuclear genome abundance across the brain, independent of the mtDNA, showed that cellular density varies by up to ~8.5-fold between brain regions (highest: cerebellum, lowest: visual cortex). Unlike in peripheral tissues, this remarkably large variation in neuronal and glial cell somata density significantly skews mtDNAcn estimates. Extended Data Fig. 2 shows the inter-correlations of mtDNAcn and mtDNA density in relation to enzyme activities. Compared to mtDNAcn, which is confounded by the presence or absence of somata/nuclear genome, mtDNA density was more consistently associated with RC enzymatic mitochondrial phenotypes across all brain regions, and therefore likely represents a more generalizable estimate of mitochondrial genome abundance across mouse brain regions.
We subsequently combined these five primary mitochondrial measures into a mitochondrial health index (MHI) by dividing RC activities (CI+CII+CIV) by mitochondrial content (CS+mtDNA density), thus creating an index of energy production capacity on a permitochondrion basis 30 (see Methods for details). Although each of the six resulting features are partially correlated with other features, the proportion of shared variance between individual features across brain regions is 31-64%, indicating that they each contribute some nonredundant information about mitochondrial phenotypes, which can be deployed at scale.
Our study design first aimed to profile animal-to-animal differences in mitochondrial phenotypes across a broad set of brain regions known to be associated with anxiety-related behavior, social behaviors, cognition, or mitochondrial disorders, and to relate these measures to each animal's behavioral phenotypes. In total, we enzymatically and molecularly phenotyped 571 samples covering 17 cortical, sub-cortical and brainstem brain regions isolated by bilateral punches at defined stereotaxic coordinates ( Table 1, Extended Data Fig. 3). To eventually compare the specificity of our findings related to brain and behavior, 5 peripheral tissues were also collected and analyzed from each animal, which we expected to show somewhat related but more modest associations with behavioral outcomes.

Protein levels and enzymatic activity
We initially explored if RC protein abundance is a viable surrogate for mitochondrial RC activity 31 , which could theoretically allow high spatial resolution imaging of the entire brain. We focused on the cerebellum due to its well-defined cellular composition and laminar organization, where the Purkinje cell layer is flanked by molecular and granular layers 14 (Extended Data Fig.   4a). Compared to protein abundance, enzymatic activity ultimately determines mitochondrial RC function and energy production capacity and consequently should be regarded as the most representative measure of mitochondrial phenotypes. In consecutive cerebellar slices, we compared RC complex II enzymatic activity measured spectrophotometrically, to the protein abundance of a complex II subunit, SDHA (succinate dehydrogenase, subunit A), for which a validated high-affinity antibody allows its quantification by microscopy (Extended Data Fig. 4bc). Across the three cerebellar layers, enzyme activity did not correlate with protein abundance assessed by immunohistochemistry and densitometry (proportion of shared variance, r 2 =0.02-0.07), indicating that protein abundance and enzymatic activity are not equivalent (Extended Data Fig. 4d). The reasons for this finding could include the action of post-translational modifications, stoichiometry of the four SDH subunits, or the biochemical context that drive biochemical activity independent of protein content (e.g., 47 ). Therefore, we focus all downstream analyses on direct measures of mitochondrial RC enzymatic activity and mtDNA density.

Diversity of mitochondrial RC activities between animals
We next examined mitochondria-behavior associations in a cohort of mice exhibiting naturally-occurring behavioral and mitochondrial variation. Our goal was to identify robust and generalizable associations between mitochondrial phenotypes and behavior, rather than potential correlations that would exist among only subgroups of animals (naïve or stressed). Therefore, to further extend the spectrum of mito-behavior variation using well-characterized rodent stress models, subgroups of the cohort were either chronically administered corticosterone (CORT) 48 for three months, or exposed to 10 days of chronic social defeat stress (CSDS) 49 (see Extended Data Fig. 3a). To create additional diversity among groups, half of the CSDS mice were allowed to recover for two months, as some stress-induced behavioral and mitochondrial changes may change again when the stressor is removed. Based on previous work 24,50 , we also expected naturally-occurring behavioral and brain molecular phenotype differences between animals that are resilient or susceptible to CSDS (based on the social interaction test; see Methods) 51 . The naturally-occurring variation in mitochondria and behavior plus the effects of various exposures provides a strong test of robustness for our hypothesis, which was that across a diverse population of mice, mitochondrial phenotypes in specific brain regions are consistently associated with behaviors.
Across our cohort of inbred male mice with a range of exposures, a wide spectrum of mitochondrial phenotypes was observed. The average variation for all measures (4 individual RC enzymatic activities and mtDNA density) across all animals and all 17 brain regions was a C.V. of 36% (coefficient of variation = standard deviation / mean). For any given brain region, the absolute variation in mitochondrial phenotypes between mice reached up to 2.9-fold between the animal with the lowest and the animal with the highest activities. This means that for a given brain region, there are large mouse-to-mouse differences in mitochondrial content and RC activities, even among inbred naïve mice not exposed to stressors (Extended Data Fig.   5). Peripheral tissues showed an average animal-to-animal C.V. of 25%, about a third less variation than for brain regions, indicating that the brain may exhibit particularly large interindividual differences.
To verify that the interventions used to enhance naturally-occurring variation in behavior and mitochondria were effective, we first compared mitochondrial features in mice exposed to CORT and CSDS relative to naïve mice. The effects of CORT and CSDS on mitochondrial phenotypes were quantified as standardized effect sizes (Hedges g) and shown in Figure 1a-c (detailed in Extended Data Fig. 6). Both interventions altered mitochondrial RC complexes, CS enzymatic activities, mtDNA density, and MHI in a region-specific manner. Although not statistically significant, CORT-treated mice tended to have higher mitochondrial activities than non-stressed animals in ~60% of brain regions, whereas CSDS animals trended towards lower activities in ~82% of brain regions, which was statistically significant for CI, CIV, and MHI measures (Figure 1d), suggesting opposing effects of these two different stress models on brain mitochondria. This could be due to either the nature of the stressor, or their durations. The amygdala (Amyg) showed the greatest CORT-induced increase in CII enzymatic activity (+49%, p<0.05, unadjusted p value), whereas the periaqueductal gray (PAG) showed the largest CSDS-induced decrease in CI activity (-42%, p<0.05, unadjusted).
To determine if brain regions were co-regulated in their stress-induced mitochondrial recalibrations, we employed a topological data analysis (TDA)-based Mapper approach. Mapper is a variant of nonlinear dimensionality reduction methods (manifold learning) that produces a graph or network embedding of the high-dimensional data (a.k.a. shape graph) while recovering projection loss using an additional partial clustering step 52,53 . When applied to the six mitochondrial features (i.e., four enzyme activities, mtDNA density, and the MHI) separately for CORT and CSDS (both measured as delta of stress vs the naïve group average), the shape graphs revealed differences in regional mitochondrial recalibrations across the two groups (Figure 1e-f). This was quantified using a graph-theoretical measure of participation coefficient (PC) 54 , where higher values of PC indicate uniformly distributed connectivity across regions (integrated network) and lower values indicate a segregated connectivity pattern. CORT-1 0 induced mitochondrial recalibrations were relatively more region-specific or segregated, whereas CSDS caused a more coherent or integrated mitochondrial response across all brain regions (Figure 1g).
Although we were not powered to make statistical comparisons between recoverysusceptible and -resilient mice based on the social interaction test for CSDS mice, we did observe relatively large differences in effect sizes (g > 0.8) in mitochondrial recalibrations between the groups that could be explored in the future work (Extended Data Fig. 7).
Additionally, compared to non-recovered CSDS mice, animals who were allowed to recover from CSDS for 57 exhibited mitochondrial phenotype changes only marginally different from control mice, and in a limited number of brain regions (Extended Data Fig. 8). This result suggests, along with the main group differences, that the brain mitochondrial recalibrations are dynamic over time scales ranging from days to weeks.
Overall, these univariate and TDA-based results established the existence of naturallyoccurring and acquired variation in brain mitochondrial phenotypes between animals, providing a strong basis to test the existence of conserved associations with behaviors.

Diversity of anxiety and depressive like behaviors
Similar to the spectrum of mitochondrial phenotypes across mice, as expected from previous work, animals also naturally exhibited large variation in their behavioral phenotypes.
Behavioral tests included potential threat ("anxiety") monitored by the open field test (OFT) and elevated plus maze (EPM), hyponeophagia monitored by novelty suppressed feeding test (NSF), and approach-avoidance conflict using the social interaction test (SI). Because specific behavioral tests are generally administered in conjunction with specific interventions (CORT, CSDS), some animals were only tested on some and not all behaviors. Both CORT and CSDS produced the expected elevation in anxiety-related behavior compared to naïve mice (Extended 1 1 (EPM and OFT, r=0.60), were not correlated (OFT and NSF, r=-0.02), or were negatively correlated (OFT and SI, r=-0.46). This indicated that each test captures different aspects of behavior, thereby providing a basis to examine how different aspects of behavior might relate to region-specific brain mitochondrial phenotypes.

Brain MHI correlates with specific behaviors
We next evaluated the extent to which mitochondrial phenotypes in different brain regions was associated with behavior across all animals (Figure 2a). Behavioral scores were transformed so that higher scores on each test indicate higher avoidance/anxiety-like behaviors To better understand which brain regions were driving the direction and magnitude in the distributions, we then examined the patterns of correlations for all 17 brain regions, for all 6 mitochondrial features, across the 4 behavioral tests (Figure 2c). In the majority of brain regions, higher mitochondrial metrics were correlated with higher avoidance scores based on OFT (average r=0.12) and EPM (average r=34), although these overall correlations did not reach statistical significance. The strongest correlations were of 0.51 for CII in the primary motor cortex (M1) and OFT (p=0.025, linear regression, unadjusted), and 0.92 for MHI in the nucleus accumbens (NAc) and EPM (p=0.0005, unadjusted). Previous research in rodents showed that mitochondrial RC function in brain regions such as the NAc is linked to complex behaviors such as social dominance and anxiety 4,31 . For NSF there was more heterogeneity between brain regions, with a similar number of positive and negative correlations. Finally, for SI, we observed the opposite relationship; animals with higher mitochondrial content and RC activities generally displayed greater sociability (lower social avoidance) (average r=-0.21, strongest for MHI in the substantia nigra (SN) and SI, r=-0.78, p=0.0035, unadjusted). Interestingly, of the six mitochondrial features, the strongest mito-behavior correlations (for 3 out of 4 behaviors) were for MHI, which may reflect the superiority of MHI as an integrative measure of mitochondrial energy production capacity over individual enzymatic and molecular features, as previously observed 30,55,56 .
As expected, the average correlation between RC activities and behaviors was significantly more consistent for brain mitochondria than for mitochondria in peripheral tissues (Figure 2d). For example, whereas mitochondrial phenotypes in several brain regions correlate with anxiety-related behavior on the EPM, mitochondrial measures in the muscles, heart, liver, or adrenal glands of the same animals on average, did not correlate with behavior. This finding aligns with recent work in mice showing that mitochondrial phenotypes exhibit strong segregation between different cell types and tissues 57 , and has two implications. First, it reinforces the specificity of these mito-behavior findings for the brain. Second, it implies that mitochondria across the brain and other tissues within an individual mouse are not equivalent, and likely differentially regulated 57 . This naturally raised the question whether specific brain regions within an animal could also exhibit independently regulated mitochondrial properties, and whether brain regions could be functionally organized into separate networks based on their mitochondrial properties.

Mitochondrial phenotype-based organization of the brain
To address this question, we first asked whether mitochondrial features in each brain region/tissues were statistically independent or correlated with other tissues. Using the same rationale that underlie functional connectivity analysis in fMRI data [58][59][60] , we first generated a correlation-based similarity matrix of all mitochondrial measures across brain regions and nonbrain tissues (Figure 3a). Within the 17 brain regions, mitochondrial features were generally positively correlated (average r=0.22, p<1 -100 , two-sample t-test), with a few exceptions. Thus, within the brain, higher mitochondrial activities in one region generally also implies higher activities in other regions (p<1 -82 , one-sample t-test).
Consistent with the co-regulation of RC enzymes within the mitochondrion, a modular structure was also apparent, indicating that mitochondrial features within each region were more similar than with other regions (p<0.0001, permutation test). Figure 3b shows the average correlations of each mitochondrial feature with other features among all brain regions -RC enzymes were most strongly co-regulated (average rs=∼0.7-0.8), followed by mtDNA density (r= ∼ 0.6), and MHI (r=∼0.5), again suggesting that MHI captures different information than individual measures. To compute how similar each region is to other regions, we computed nodal degree (i.e., average inter-regional correlation, similar to the concept of global connectivity) for all regions. Nodal degree was highest in the cerebellum (average r=0.34) and lowest in the brainstem vestibular nucleus (r=0.05) (Figure 3c), suggesting that the region-to-region similarity in mitochondrial phenotypes is not randomly distributed across the mouse brain.
As suggested above from the divergence between brain regions and peripheral tissues in the mito-behavior correlation patterns, the mitochondrial features in peripheral tissues were not correlated with brain mitochondria (average r=0.02) or with other peripheral tissues mitochondria (average r=-0.03) (Figure 3d). The lack of association between brain and nonbrain tissues of the same animals suggest that MHI (as well as content and RC activities, data not shown) is not a ubiquitous animal-level feature, but rather relatively independently defined in the brain, and further specified among each brain region.
Given the overall positive connectivity across the brain and recent evidence of circuitlevel metabolic coupling across large neural networks in the Drosophila brain 12 , we examined whether certain regions exhibited particularly strong mitochondrial feature co-regulation. Functional brain networks in the living brain are defined by synchronous activity patterns (e.g., To examine whether this functional organization of the brain revealed using mitochondrial phenotypes was also evident cross-modally, we used two Allen Brain Atlas datasets of brain-wide gene co-expression 62 , and EYFP-labeled axonal projections that define the structural connectome 63 , developed in the same mouse strain. Specifically, we examined whether gene expression correlations and structural connectivity within brain regions that are functionally grouped together as a network, based on mitochondrial features, is higher than expected by chance. We used two independent graph theoretical metrics: strength fraction (S.F.) 64 and quality of modularity (Q_mod) 65 , to examine whether similar communities exist in gene and structural connectivity data as in our mitochondrial phenotype data. Using permutation testing, we randomly shuffled community structure of brain regions 10,000 times to determine whether the community structure derived from mitochondrial features is also evident in gene coexpression and structural connectivity data. Both strength fraction and quality modularity statistics indicated that mitochondria-derived networks also have higher similarity in gene coexpression (S.F. p=0.020; Q_mod p=0.008) and structural connectome data (S.F. p=0.029; Q_mod p=0.015) than expected by chance (Extended Data Fig. 11). Hence, this provides convergent multimodal evidence of mitochondrial phenotypic organization overlapping with gene expression and structural connectivity.

Network-level mitochondria-behavior correlation
Finally, to examine the potential significance and added value of this effective brain-wide mitochondrial connectivity in relation to animal behavior, we used the mitochondria-derived network results to partition the brain into three networks, and then analyzed the ability of each network to linearly predict behaviors. The mitochondrial-behavior correlations varied both in strength and direction between the three networks, and the mitochondrial phenotypes among the three networks had largely divergent associations with behaviors. The average mitochondrial activity in the cortico-striatal network consistently showed the strongest average correlations with behaviors measured in the OFT, EPM, and SI tests (r=0.20, 0.54, -0.41, respectively). For the EPM, mitochondrial features in the cortico-stratal network (1) accounted for up to 33% of the animal-to-animal variance in behavior (EPM). In contrast, the salience/spatial navigation network (2) exhibited the weakest average correlations for the same three behaviors (r=-0.03, 0.18, -0.14), accounting for <5% of behavioral variance (Figure 3g).
Together, these results suggest the presence, specificity, and functional significance of brainwide mitochondrial networks, embedded within existing neural circuitry in the mouse brain.

Discussion
Using a high-throughput approach to functionally phenotype hundreds of brain samples from a heterogenous mouse cohort, we have defined brain-wide associations between mitochondrial phenotypes and behaviors. Combined with previous findings 31, 35 , the diverging mito-behavior associations between brain regions, and between brain and non-brain tissues, brought to light the possibility that different brain regions might exhibit different mitochondrial phenotypes. In particular, the network characteristics of mitochondrial phenotypes across the brain provided evidence for the modular distribution of mitochondria across cortical and subcortical regions, as well as their relevance to animal behaviors. Based on these data, we conclude that mouse brain mitochondria may exist as behaviorally-relevant networks overlapping with, but distinct from, other modalities including gene expression and structural connectivity.
Utilizing CORT and CSDS interventions to induce behavioral variation in our animal cohort also allowed us to compare the effects of these two interventions on mitochondrial phenotypes. Some brain regions were found to respond, in some cases, in opposite directions, particularly after exposure to CORT. In contrast to CORT, the recalibrations of brain mitochondria to CSDS was more uniform, with the majority of brain regions exhibiting a coordinated reduction in most mitochondrial features. This difference in mitochondrial recalibrations between both stress models may be driven by several factors. This includes the stressor duration, although the temporal dynamics over which stress-induced mitochondrial recalibrations take place remain porrly defined and will require further focused attention.
Differences in the effects of CORT vs CSDS on mitochondria also could be related to their neuroendocrine underpinnings (single hormone for CORT vs multiple physiological neuroendocrine recalibrations for CSDS), regional differences in glucocorticoid and mineralocorticoid receptor density, or other factors generally relevant to interpreting chronic stress rodent models.
One valuable aspect of our study is that it highlights the specificity of earlier findings indicating a connection between mitochondrial RC function in the NAc and anxiety behaviors 4,31 .
In our study, the strongest correlation between EPM-based anxiety-like behavior and MHI was in the NAc, confirming the strong association between NAc mitochondrial energy production capacity and anxious behavior. We also extend this finding to show that mitochondrial phenotypes are linked to behaviors across not only isolated regions, but likely distributed among brain networks. Similar to the conceptual shift from regional and cellular perspectives towards distributed brain networks, circuits, and neuronal ensembles 66,67 , our findings therefore advance the notion that mitochondria may modulate brain function and behavior through distributed mitonetworks. This notion is consistent with recent evidence of metabolic coupling with distributed brain-wide patterns of neural activity linked to behavior in Drosophila 12 .
While we cannot directly explain why distinct mitochondrial phenotypes appear to exist across brain regions, this may be driven by three main factors. First, mitochondria could respond to differences in neuronal circuit functioning, in agreement with observed coupling of neuronal and metabolic activities 12 , such that the cellular infrastructure of neural circuits that fire together not only wire together, but also generate similar mitochondrial phenotypes (i.e., mitotypes). A second possibility is that brain regions which are regularly co-activated, i.e., within the same functional networks, harbor similar levels of receptors for neuroendocrine factors (stress and sex steroids) known to influence mitochondrial biogenesis and/or functional specialization 68,69 . A third possibility would involve differences in metabolites and substrates that regulate mitochondrial RC activities, arising similarly among co-activated brain regions. These and other potential factors underlying the modularity of mitochondria within the mouse brain remain to be investigated. showing that quantifying mtDNA content on a per-cell basis (mtDNA:nDNA ratio) is heavily skewed by cellularity variations across brain regions and not directly related to RC energy production capacity. As such, the mtDNA:nDNA ratio (mtDNAcn) is driven by how many cell bodies are present in the tissue, and correlates poorly with either mitochondrial content or RC activity. Therefore, our data reinforce the notion that mtDNAcn on its own is not a valid measure of mitochondrial phenotypes 73,74 . In the mouse brain, our findings suggest that mtDNA density per unit of volume (rather than per cell) is a more biologically meaningful mitochondrial feature when comparing brain regions that differ in cellularity.
By examining mitochondrial features across the mouse brain, we discovered a moderate level of global functional connectivity across mitochondria in most brain regions, but not among peripheral tissues (Figure 3). By functional connectivity, we do not imply that mitochondria are directly connected to each other in the same way that neurons project and chemically (de)activate each other; but rather that they share functional properties. If mitochondrial phenotypes are directly determined by genetic and physiological factors, the logical expectation is that all mitochondria within different organs of the same organism should exhibit a high degree of coherence (i.e., correlated with each other). In other words, the animal with the highest mitochondrial content or RC activities in the one brain region should also be the animal with the highest activities in other brain regions and tissues. Our results strongly disprove this point. Similarly, previous work on multiple human tissues showed that mtDNAcn was not significantly correlated between organs, including across three brain regions 75 . Here we extend these data to 17 brain regions and 5 non-brain tissues, demonstrating that identified correlations are relatively modest (the strongest is r=0.31 for the cerebellum, representing less than 10% of shared variance). Notably, while on average we observed some degree of mitochondrial phenotype coregulation across the brain, some brain regions exhibited no consistent correlation with other regions. The most parsimonious explanation for this result is that individual animals differentially recruit different circuitry, which secondarily shape their mitochondria and drive animal-specific patterns of regional mitochondrial variation. Considering theories that stipulate that the brain strives for maximal energetic efficiency 1 and that there may be limits to cellular energy conversion rates 76 , we also cannot exclude the possibility that a limited quantity of resources (i.e., mitochondrial content or energy conversion) is available within the brain or the whole organism, which are distributed unequally in an activity-dependent manner among different organ systems, functional networks, and individual brain regions. Thus, we speculate that different animals may achieve an optimal balance of systemic and neural functions through specific combinations of mitochondrial activity in different brain regions, a possibility that remains to be tested.
In trying to better understand why certain regions displayed stronger brain-wide connectivity than others, we identified mitochondria-based communities, or networks of brain Until now, the methods used to reveal neurobiological and metabolic networks in mammals have typically been through indirect functional and structural connectivity analysis 20 .
Here we have developed a scalable approach to examine mitochondrial phenotypes across a large number of brain regions in mice with a range of behavioral phenotypes. We showed that mitochondrial phenotype connectivity was non-random, and linked with gene co-expression and structural connectivity, thereby providing converging evidence of mitochondria-based networks across modalities. This study synergizes with recent work 12 providing the technical and empirical foundation to bring mitochondrial biology into brain-wide, network-based models of neural systems in mammals. Applied to mitochondria, network-based analytics should contribute to develop more accurate maps, and eventually an understanding of what drives mitochondrial and metabolic properties across interacting neural circuits. Developing a spatiallyresolved understanding of brain mitochondrial biology will help to resolve the energetic contraints on brain function and behavior.
Limitations. Notable limitations of this study include the lack of cell type specificity.
Neurons operate in a metabolic partnerwhip with astrocytes and glial cells 81 , and different cell types exhibit different molecular mitochondrial phenotypes (e.g., 82 ) that cannot be disentangled in tissue homogenates. While neither enzymatic or functional mitochondrial profiling at the single brain cell is currently technically feasible, it remains possible that mitochondrial phenotypes between brain regions are influenced by differences in cell type proportions.
Moreover, our molecular and biochemical mitochondrial phenotypes do not reflect other factors that can influence the efficiency or activity of the mitochondrial RC or OxPhos system in vivo, such as variations in cofactor abundance and RC structural assembly (e.g., supercomplexes) 83 .
Because functional assays require harvesting brain tissue, it also was not feasible to asertain within a given animal how stable (trait) or dynamic (state) the brain biochemical mitochondrial phenotypes are. If mitochondrial phenotypes were more dynamic than expected, our estimated proportion of behavioral variance attributable to mitochondrial biology could be substantially underestimated.

Animals
This study was carried out in accordance with NIH Guidelines, and was approved by the Institutional Animal Care and Use Committee (IACUC) at New York State Psychiatric Institute.
Three month-old CD1 retired breeder male mice were obtained from Charles River, and were used as the aggressors in the social defeat model.

Chronic Social Defeat Stress (CSDS)
Aggressors Prescreening: All CD1 mice used in the experiment were pre-screened for aggressive behaviors as previously described 50 . During a three-day screening procedure, a novel C57BL/6J mouse was placed in the cage of the CD1 mouse for 3 min. C57BL/6J screener mice were not further used in the study. The latency of the CD1 mouse to attack the C57BL/6J screener mouse was recorded. CD1 mice that attacked in less than 1 minute on at least the last two consecutive screening days were considered to be aggressive.
Experimental Groups: Thirteen C57BL/6J mice underwent social defeat stress, while 6 remained in their cages, to serve as the naïve group. Animals were randomly assigned to these groups. One mouse subjected to social defeat died of unknown causes during the duration of the experiment. Of the remaining 12 social defeat mice, half (n=6) were sacrificed two days after the completion of the stressor, and they are referred to as the 'stressed' CSDS group. Three naïve mice were sacrificed at the same time. The other half of the CSDS mice (n=6) along with 3 naïve mice, were allowed an 8.5 weeks (59 days) stress recovery period prior to sacrifice.
This group is referred to as the 'recovered' group. Animals were selected for the recovery group from the CSDS group so that there was an even ratio of susceptible to resilient mice in the 'stressed' and 'recovered' groups.  Behavioral Z-Scores: Social avoidance scores were determined by averaging the z scores for 4 measures as previously described 24 ; SI ratio, corner ratio, time spent in SI zone, and time spent in corner zones. The z scores for SI ratio and time spent in SI zone were multiplied by -1, so that higher scores across all 4 measures indicated higher avoidance.
Similarly, the OFT score was determined by averaging the z scores for the 2 measures. The z scores were multiplied by -1, so that higher OFT score indicate higher avoidance of the brightly lit center of the OF. EPM scores were also multiplied by -1, so that higher EPM scores indicate higher avoidance of the open arms of the EPM. High NSF scores already indicate higher avoidance, so they were not inverted. Therefore, across all 4 tests, higher scores indicate higher anxiety-like/avoidant behavior.

Tissue Collection
Animals were sacrificed by rapid decapitation to maintain mitochondrial integrity in brain tissue. Brains were rapidly flash frozen in ice-cold isopentane, stored at -80°C, and later transferred to -170°C (liquid nitrogen vapor) until mitochondrial measures were performed.
Brains were transferred back to -80°C during the week prior to being sectioned and were then transferred to -30°C the night before sectioning. Brains were sectioned coronally on a Leica Model CM3050 S cryostat. The cryostat internal temperature and blade temperature were set to -22°C during sectioning. Brains were mounted onto a specimen disk using optimal cutting temperature (OCT) compound. The brain was sectioned coronally, alternating between two 200μm thick slices and then two 20μm, and were deposited onto microscope slides. Brain sections were kept at -80°C until brain region-specific tissue sample collection. One of the recovered CDSD mouse brain contained blood and could not be reliably sectioned, so this animal's brain tissue was excluded from analysis. A second CSDS mouse brain and a naïve brain both cracked during slicing, so some of the brain regions (n= 7 for CSDS mouse, n=8 for naïve) could not be obtained from those brains, but the regions that were obtained were included in analysis.
Tissue biopsy punches on frozen brain sections: The scalable Allen Mouse Brain volumetric atlas 2012 86 was used to determine the location of each brain region of interest, and their distance from bregma. The atlas' Nissl-stained images were used as a reference for estimating each section's distance from bregma, which determined which brain slices would be used for punch biopsy collection of each region of interest (Extended Data Fig. 3b). Brain slices were punched using 1.00mm diameter Robbins True-Cut Disposable Biopsy. Two bilateral punches were collected for each brain region over dry ice. For brain regions in the midline, punches were taken from two consecutive slices. All tissue punch locations were approximated by using landmarks on the slice and by comparing to the atlas. It is important to note that the 1mm puncher was larger than the actual brain region in some instances, and so parts of neighboring brain regions may have been included in the punches. Therefore, the punching technique may include some error. Punches were stored at -80°C until they were ready to be used for enzymatic activity assays. Intact consecutive 200μm cerebellar slices were stored at -80°C for immunohistochemical staining.
The tissue punches were too light to be accurately weighed, so we had to estimate the weight based on the reported mouse brain density of 1.04g/cm 387  we obtain an estimated mass of 0.163mg per punch, thus 2 punches were approximated to weigh 0.327mg.

Mitochondrial Measurements
Tissue Preparation: The punches from each brain region were homogenized in 0.2mL of homogenization buffer (1mM EDTA and 50mM Triethanolamine) (2x 1mm tissue punch/ 0.2mL of homogenization buffer), with 2 tungsten beads to disrupt the tissues' cells and release the mitochondria. Tissues were homogenized using a Tissue Lyser (Qiagen cat# 85300), which was run at 30 cycles/sec for 1 minute. The tissues were then incubated in ice for 5 minutes, and were then re-homogenized for 1 minute. Tissues were vortexed to ensure homogeneity.
Peripheral tissues were cut over dry ice, weighed, and were then homogenized 1:180 (weight:volume, mg: μ L), except for heart samples that were further diluted to 1:720 to be in the dynamic range of the assays.

Enzymatic activities
Enzymatic activities were quantified spectrophotometrically for Citrate Synthase (CS), complex I (CI, NADH-ubiquinone oxidoreductase), Succinate Dehydrogenase (CII, succinateubiquinone oxidoreductase, also known as SDH), Cytochrome C Oxidase (CIV, COX) and were expressed per mg of tissue, as described previously 32 , with some modifications as described below in full details. All miniaturized assay measurements were performed in 96-well plates and enzymatic activity assays recorded on a Spectramax M2 (Spectramax Pro 6, Molecular Devices). Linear slopes reflecting changes in absorbance of the reporter dye were exported to Microsoft excel and converted into enzymatic activities using the molar extinction coefficient and dilution factor for each assay. The assays were optimized to determine the minimal amount of tissue required to obtain reliable results assessed by the C.V. between duplicates. The assays were then further optimized to determine the minimal amount of brain tissue homogenate required for each individual assay. Assay validation involved regressing increasing tissue amounts (number of punches) with observed activities, which confirmed that an increase or decrease in tissue used produced proportional changes in total activity.
To validate the miniaturized biochemical enzymatic assays, we performed each assay in both the miniaturized 200ul format in 96-well plates, and the traditional 1ml cuvette format, using increasing tissue homogenate volumes (4, 8, 12, 16, 20ul) from the same brain region (cerebellum) of a wild-type control mouse. The same homogenized tissue (1:200 weight:volume, mg:μL) was used for CS, CI, CII, and CIV spectrophotometric assays from which the respective enzyme activities were quantified using the reagents and procedures described below. Both the miniaturized and standard-size cuvette assays showed high agreement (r 2 = 0.81-0.96; ps = 0.0032-0.037) (Extended data Fig. 1 Samples were run in duplicates for each enzyme, along with a nonspecific activity control, and every plate had a positive control (heart homogenate). The 96-well plates were designed so that each brain region/tissue from all animals were run on a single plate, which prevents potential batch variation for comparisons between the animals. However, due to the size of the plates, no more than 2 types of tissues could be run together on a single plate. This

Confocal imaging
The immunofluorescent staining was imaged on a Leica SP8 confocal microscope

SDHA staining analysis
Final images were deconvolved in the Lightning software using the default settings.
Images were analyzed in ImageJ version 1.52p. Fifty consecutive slices from each stack were selected for a section thickness of 7.5μm, which roughly covered the depth of penetration of the antibodies. The 16bit images were first thresholded using the Triangle algorithm within the integrated Auto-thresholding plugin to create a binary image. The particle analyzer was then used to quantify the percent of the area stained for each slice. The percent area stained from each slice was summated to calculate the percent of volume stained. For compartmentalized data, the molecular and granular layers were separated using a rectangular selection of approximately 50px width through each stack. For the purkinje layer, a rectangle was drawn that narrowly captured the stained area of each slice. The volumetric staining for each region of interest was calculated from the percent of area stained quantified using the particle analyzer.

TDA-based mapper analysis
After creating a delta of mitochondrial features between stressed (CORT or CSDS) and naïve mice, the data matrix for each group was processed through the TDA-based Mapper pipeline 52 . The input data matrix contained 102 concatenated rows for brain regions (17 brain For example, linear filter functions like classical principal component analysis (PCA) could be used to preserve the global variance of the data points in the high dimensional space. However, a large number of studies using animal models and computational research suggest that interregional interactions in the brain are multivariate and nonlinear [89][90][91] . Thus, to better capture the intrinsic geometry of the data, a nonlinear filter function based on neighborhood embedding was used 52 . The second step of Mapper performs overlapping n-dimensional binning to allow for compression and reducing the effect of noisy data points. Here, we divided the data into lower dimensional space into 64 bins with 70% overlap. Similar results were observed for different number of bins (e.g., 49, and 81). Third, partial clustering within each bin is performed, where the original high dimensional information is used for coalescing (or separating) data points into nodes in the low-dimensional space. Partial clustering allows to recover the loss of information incurred due to dimensional reduction in step one 92 . Lastly, to generate a graphical representation of the shape of the data, nodes from different bins are connected if any data points are shared between them.
The Mapper-generated graphs can be annotated (or colored) using meta-information that was not used to construct the graphs. Here, we annotated these graphs using region-labels to examine whether mitochondrial features were similarity expressed across all regions or whether regional specificity was observed across the two groups. To quantify the extent of segregation (high regional specificity) or integration (low regional specificity) across the six mitochondrial features, we used a graph theoretical measurement of participation coefficient 54 .
Participation coefficient ܲ of a node ݅ is defined as: is the number of links of node uniformly distributed among all communities of the graph (and hence integrated) and it is close to 0 if its links are mostly within its own community (and hence segregated).

Multi-slice community detection
One of the most commonly studied mesoscale aspect of a graph is modularity, where highly modular graphs consist of cohesive groups of nodes (or communities) that are more strongly connected to each other than they are to the rest of the network 93 . Examination of modularity has been recently extended to multi-slice networks that are defined by coupling multiple adjacency matrices across time or modality 61 . Here, we used six slices, each derived from one of the mitochondrial features, where each slice represented weighted adjacency matrix of pair-wise Pearson's correlation between brain regions. We use categorical multi-slice community detection algorithm with the presence of all-to-all identity arcs between slices 61 . The generalized modularity function for multi-slice community detection is given by represents weighted adjacencies between nodes (1,000), followed by consensus clustering to get stable results for identifying large-scale networks of brain regions 97 .

Comparison of mito-based networks with transcriptomic-and structural connectivitybased networks
To compare our mitochondrial communities against other modalities, we examined the organization of gene co-expression in the mouse brain as well as mouse structural connectome data. We utilized the Allen Brain Atlas' ISH (in situ hybridization) feature, which maps each gene onto a reference standard coordinate atlas image, providing a spatial estimate of transcript levels representing gene expression 62 . Specifically, we used the Anatomic Gene Expression Atlas (AGEA) (https://mouse.brain-map.org/agea), which is a a three-dimensional male adult C57BL/6J mouse brain atlas based on the ISH gene expression images. The AGEA feature allows for selecting both a 'seed voxel' and a target 'selected voxel' from exact brain coordinates coordinates, yielding the transcriptome-wide correlation between the seed and target regions. The correlation is a measure of the average co-expression between the two voxels. The co-expression values were obtained for all possible pair of brain regions (17x17 matrix), yielding a gene co-expression correlation matrix to which the structure of our mitobased networks could be compared.
The structural connectome data were obtained from The Allen Mouse Brain Connectivity Atlas 63 , which is a mesoscale connectome of the adult mouse brain. We utilized the normalized projection strength values for the all brain regions of interest. Because this atlas does not distinguish between dorsal and ventral dentate gyrus, data was obtained for 16 brain regions.
For three brain regions where the Atlas provides connectivity values for multiple subregions (OFC, VN, Cerebellum), the connectivity values were averaged across the subregions to yield a global measure for each region, thus matching the dimensionality of our mitochondrial dataset.
We used non-parametric permutation statistics with 10,000 permutations to examine whether the mitochondrial features derived networks were also more densely connected than expected by chance in other modalities (transcriptomic and structural connectivity data). To measure the degree of within-network connectedness we used two established metrics: modularity index (Q_mod; 65 ) and strength fraction (S.F.; 64 ). As presented in Extended Data  Fig. 11, for each modality, and across the two metrics, the networks derived from mitochondrial features were more tightly knit than expected by chance (all ps<0.05). Thus providing a convergent multimodal evidence of mitochondrial phenotypic organization overlapping with gene expression and the structural connectome.

Statistical analysis
Standardized effect sizes (Hedge's g) were computed to quantify the effect of stress conditions on mitochondrial features. Significant effect sizes were determined by the 95% confidence intervals. Two-way ANOVAs were used to compare the effects on a variable between groups. Frequency distributions of the effect sizes for the two stressors and of the mitochondrial-behavior correlations were fitted with Gaussian curves and analyzed by one sample t tests compared to the null distribution. A survival curve using Mantel-Cox log-rank test was generated for novelty suppressed feeding test latencies because the test is capped at 600 sec. Behavioral scores between groups were analyzed using Tukey's multiple comparison ordinary one-way ANOVA. Correlations between behavioral scores and mitochondrial activities were estimated using Spearman's r (r) to guard against inflation. Correlations between tissues' mitochondrial activities were measured as Pearson's r, and hierarchical clustering was performed on the data using Euclidian distance with Ward's clustering algorithm. T-tests were used to compare groups' correlations. Permutation testing was used to analyze mitochondrial features between versus within brain regions and to analyze the mitochondrial networks against gene co-expression and structural connectome data. To assess the group (CORT vs CSDS) differences in the regional response to stressors are statistically robust, a phase randomized (PR) null approach (Theiler et al. 1992) was used, where we generated 1000 iterations of null data separately for CORT and CSDS groups and ran our TDA-based Mapper approach on each null dataset. PR null preserved the covariance across mice but shuffled any relation between regional responses. Non-parametric permutation statistics were later estimated for each group.

Data Availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author upon request.

Code Availability
For   Nodes in the Mapper output represent regional mitochondrial features (rows from input matrix) that are highly similar across mice. Thus, the brain regions that undergo similar stress-induced recalibrations in specific mitochondrial features cluster together in single nodes (most similar) or interconnected nodes (moderately similar), whereas regions that undergo divergent recalibrations are not connected. The pie-chart-based annotation of graph nodes allowed us to examine the degree of co-regulation of mito-features across brain regions. (f) Topological data analysis (TDA) based Mapper approach to determine if brain regions were co-regulated in their stress-induced mitochondrial recalibrations for the two groups. (g) Graph-theoretical measure of participation coefficient (PC), measuring the degree of segregation/integration among brain regions in each graph, indicating that CORT-induced mitochondrial recalibrations were more region-specific (or segregated) whereas CSDS caused a more integrated response. A phase randomized null approach was used to estimate the robustness of these results, such that null data were generated where covariance across mice (columns in the input matrix) was preserved while the regional relations (rows) were disrupted. Non-parametric permutation statistics revealed CORT group has a significantly lower participation coefficient (PC) as compared to the Mapper graphs generated from the null models (p=0.037). A similar analysis for the CSDS group revealed a trend towards significance for higher PC as compared to the null models (p=0.069).  Fig. 8). (c) Individual correlations for the 17 brain region, across the 6 mitochondrial features, for each behavioral score, quantified as Spearman's r. OFT and EPM behavioral scores were inverted so that higher scores on all four tests indicate higher anxiety (see Extended Data Fig. 8 for additional details). The strongest correlations for each behavioral test are denoted by yellow boxes, with the scatterplots shown below. An adjusted p value of <0.002 was applied (false-discovery rate 1%). All tests have been adjusted so that a higher score indicates higher anxiety-like behavior (see Methods for details). (d) Average correlation of each behavior for the brain (B) and tissue (T) mitochondrial features; two-way ANOVA. n=10-27 mice per behavioral test. Correlation between the average mtDNAcn and average mtDNA density per brain region. (c) Correlations between mtDNAcn (left) and mtDNA density (right) with the other five mitochondrial measures in all animals for each brain region and tissue, measured by Spearman's r. The strongest positive and negative brain correlations for each are highlighted with yellow boxes and plotted in (d). Fig. 3 Fig. 7. Exploratory analysis comparing susceptible and resilient subgroups of CSDS mice. Effect sizes for mitochondrial outcomes shown separately by brain region for animals classified as susceptibe or resilient to CSDS (Methods). Effect sizes are Hedge's g, with significant effect sizes (95% confidence interval) labeled with the % change compared to naïve (non-stressed) mice. Fig. 8. Effect of recovery from social defeat stress on mitochondrial features across brain regions. (a) Average effect of 2 months recovery (Rec) from CSDS stress across the 6 mitochondrial features for each brain region, as compared to non-recovered CSDS mice. Effects sizes are quantified by Hedges g, with 95% confidence intervals. MHI in recovered mice was higher than non-recovered CSDS in some regions (Caudoputamen, Hypothalamus-CA3), and lower in others (Hypothalamus, Medial orbitofrontal cortex, Vestibular 1 nucleus). (b) Average effect sizes of 6 mitochondrial features of non-recovered CSDS mice and recovered mice compared to naïve mice in three selected brain regions from A (red boxes). (c) Frequency distribution of the effect size of CORT, CSDS, and recovery as compared to naïve mice on all 6 mitochondrial measures in all 17 brain regions, gaussian-fitted curve, indicating that recovered mice have lower MHI than naïve mice and marginally higher MHI than CSDS (****p<0.0001, one-sample t test (two tailed)). (d) TDA analysis using Mapper of mitochondrial measures in naïve mice, CORT, CSDS, and CSDS-recovered mice, revealing the recovered group to be more similar to the non-recovered CSDS group than to naïve mice. Fig. 9. Behavioral test results (a) To account for both the % time and % distance traveled in the center, the two measures were z-scored and averaged, creating a single Open Field Test (OFT) score. Z-scores were inverted so that a higher OFT score represents less time and distance in the center, indicating higher anxiety. Open field was run for all groups. (b) Elevated plus maze (EPM), measuring percent of time spent in open arms/total time, with the score being multiplied by -1 so that higher scores represent less time spent in the open arms, which indicates higher anxiety. EPM was run only on CORT mice. (c) Novelty suppressed feeding (NSF) test used to measure the latency to feed in a novel environment, with the test being capped at 600 sec. Cum. summ=cumulative survival, meaning percent of mice that have not eaten; Survival curve: Mantel-Cox log-rank test. NSF was run only on CORT mice. (d) Social Interaction test (SI), represented by a Social Avoidance score for all CSDS mice, separated by susceptible (sus.) and resilient (res.), and recovered (rec). Social avoidance was measured as the average of 4 scores; z scored social interaction ratio, z scored time spent in interaction zone, z scored time spent in corner, and z scored corners ratio. Social interaction ratio and time spent in SI zone were inverted so that higher values for all four measures indicate higher avoidance. p<0.05, p<0.01, Tukey's multiple comparison ordinary one-way ANOVA. Fig. 10. Frequency distribution of all correlations of mitochondrialbehavior pairs by behavioral test. Gaussian fits, p values from one-sample t test (two tailed) againt the null hypothesis (r=0). Fig. 11. To examine whether the functional organization of the brain revealed using mitochondrial features is also evident cross-modally, we utilized non-parameteric permutation stastics with 10,000 permutations to compare the topology of our mito-derived networks with the expected topology based on the transcriptomic and structural connectomes in the Allen Brain Atlas (see Online Methods for details). Specifically, we examined i) whether the tightly knit networks (or communities) observed in the mitochondrial features are more densely connected than expected by chance (a-b); and ii) whether the networks derived from mitochondrial features were also more densely connected than expected by chance in other modalities, including gene-coexpression data (c-d) and EYFP-based structural connectivity (ef). To measure the degree of within-network connectedness we used two established metrics: modularity index (Q_mod; Newman 2006) and strength fraction (S.F.; Richardi et al. 2015). The histograms depict distribution of within-network connectedness generated using 10,000 permutations (or data shuffling). The real value of within-network connectedness is shown using red dots across all plots. For each modality, and across the two metrics, the networks derived   Mito feature specifically recalibrated in brain Region C   Average Inter-regional Correlation

Global connectivity by brain region
Average inter-regional correlation Average inter-tissue correlation Chance level Regional co-variation for each mito feature across 17 brain regions 6 mito features

Mitochondria-driven community detection among brain regions
Categorical multi-slice community detection